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Synopsis 

The  paper  presents  an  assimilation  of  mathematical  models  and  solutions 
needed  in  order  to  develop  computer  based  analysis  of  dynamic  structures. 
Using  the  variational  formulation  and  a  direct  integration  technique,  a  dynamic 
finite  element  model  is  developed.  Modal  analysis  of  unknown  displacements  of 
the  structure,  and  the  dynamic  reduction  of  the  structure  are  presented  as 
alternative  solutions.  A  system  of  micro-computer  based  programs  which  apply 
the  presented  solution  techniques  is  described.  The  system  of  programs  support 
varying  cross  sections  of  frame  members,  application  of  static,  harmonic  and 
non-harmonic  loading  conditions,  and  node  displacements^  in  the  form  of 
uniform  base  motion  or  independent  node  movement.  V>  / 
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Introduction 

Computer  analysis  of  dynamic  structures  has  for  some  time  been  limited  to 
mainframe  computers.  The  importance  of  conducting  a  detailed  analysis  of  any 
structure  is  evaluated  against  access  to,  and  the  cost  of  using  a  mainframe 
application  to  do  that  analysis.  There  are  situations  where  analysis  by 
mainframe  is  not  possible  or  is  not  justified.  In  such  cases,  solution  by  hand 
may  be  impractical.  There  is  a  need  to  conduct  rigorous  analysis  of  dynamic 
structures  that  are  too  simple  to  justify  using  mainframe  applications  and  too 
complicated  to  be  solved  by  hand.  Micro-computers  are  viewed  as  a  possible 
means  of  satisfying  this  need. 

The  large  amounts  of  memory  required  by  the  techniques  which  enable  dynamic 
structures  to  be  modeled  in  a  form  solvable  by  a  digital  computer  have  restricted 
their  implementation  on  micro-computers.  However,  these  techniques  continue 
to  be  studied,  refined,  and  combined  with  other  techniques  in  the  attempt  to 
develop  an  optimal  solution.  In  addition,  micro-computers  with  abilities  to 
address  memory  measured  in  the  multi-megabytes  are  becoming  widely  available. 
With  improved  techniques  and  larger  memory  capacities,  one  can  expect  that 
rigorous  analysis  of  simple  dynamic  structures  will  soon  be  done  conveniently 
and  inexpensively  using  micro-computers. 

Toward  that  end,  the  mathematical  formulations  required  to  model  dynamic 
structures  on  a  micro-computer  are  synopsized.  Combining  these  methods  with 
a  technique  of  reducing  the  complexity  and  number  of  resulting  equations  then 
results  in  an  useful  engineering  analysis  tool. 

The  paper  first  illustrates  how  the  Finite  Element  Method  is  used  to  discretized 
the  problem  and  express  it  in  a  matrix  form.  Next,  the  Newmark  method  of  direct 
integration  is  used  to  simplify  resulting  integrations  with  respect  to  time. 
Further  simplification  of  the  equations  are  made  possible  through  formulation 
and  solution  of  the  eigenvalue  problem.  Finally,  a  method  for  reducing  the 
number  of  equations  which  must  be  solved  is  presented. 

To  show  how  the  above  techniques  are  applied  to  a  micro-computer,  a  system  of 
programs  is  described.  The  programs  are  capable  of  solving  the  resulting 
equations  for  dynamic  analysis  of  undamped  linear  plane  frame  structures  using 
any  of  the  presented  solutions.  Flow  diagrams  and  program  listings  are  provided. 
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Finite  Element  Formulation 

The  Finite  Element  method  is  widely  used  in  the  analysis  of  structures.  It  has 
the  ability  to  systematically  describe  a  structure  in  a  matrix  form  which  is 
easily  applied  to  computer  computation.  Understanding  the  methods  by  which  the 
matrix  form  is  arrived  at  is  important  in  understanding  the  capabilities  and 
limitations  of  a  computer  application  which  employs  the  method. 

The  derivations  presented  in  this  section  draw  heavily  from  a  text  by  J.  N.  Reddy, 
"An  Introduction  to  Finite  Element  Method"  (see  the  bibliography). 

Discretization 

This  model  will  describe  a  structure  as  an  assemblage  of  two  node  frame 
elements.  Each  node  will  have  three  degrees  of  freedom,  horizontal,  vertical, 
and  rotational  movement.  Mathematically,  the  frame  element  will  consist  of  a 
superimposed  one-dimensional  bar  element  and  a  two-dimensional  beam  element. 
The  bar  and  beam  element  are  superimposed  in  a  manner  that  assumes  the 
transverse  and  rotational  deflections/loads  are  independent  from  axial 
deflections/loads. 

Bar-Element 

The  governing  differential  equation  for  the  bar  element  is: 

82u  8 

m  — -  ♦  — 

8t2  8x 


8u 
AE  — 
8x 


♦  F(x,t)  =  0 


Where  F(x,t)  is  an  axial  forcing  function  which  varies  linearly  with  x,  m  is  the 
mass  per  unit  length,  A  is  the  cross  sectional  area,  and  E  is  the  modulus  of 
elasticity.  Damping  has  been  ignored. 


The  variational  formulation  is  found  by  integrating  the  governing  equation 
against  a  test  function  over  h,  the  length  of  a  bar  element. 
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Integrating: 


fh 

.  - 

32u  3v  3u 

3u 

vm  —  ♦  AE - +  vF(x,t) 

dx  -  vAE — 

3t2  dx  3x 

dx 

0 

,  « 

The  last  term  of  the  above  expression  corresponds  to  the  natural  boundary 
conditions  at  either  end  of  the  bar  element  and  will  be  denoted  as  P1f  the  axial 
force  on  the  left  side  and  P2.  the  axial  force  on  the  right  side  of  the  element. 

Assume  that  u  is  interpolated  by  a  linear  expression  of  the  form: 

2 

U=  £  Uj(t)  ^(x) 

j=l 

Assuming  that  u  and  t  can  be  separated,  for  any  given  time  t>0,  the  above 
expression  is  substituted  for  u,  and  v  =  f  j(x).  The  matrix  formulation  results: 

IMMu"}  ♦  IKMu}  =  (F(t)} 

Where  (’)  means  differentiation  with  respect  to  t  and: 


fh 

Fj  =  ♦jF(x.t)  dx  ♦  Pj(t) 

J0 

Note  that  in  the  physical  meaning  of  the  above  expressions,  Mjj  Kjj  do  not  vary 
with  time.  While  Fj  and  Pj(t)  vary  with  time,  solutions  will  be  based  on  the 
specific  values  of  F(x,t)  and  Pj(t)  at  given  points  in  time. 
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The  interpolation  functions  «J»j  (for  i=1  to  2)  must  be  sufficiently  differentiable, 

independent  of  one  another,  complete,  and  must  satisfy  the  essential  boundary 
conditions.  The  expressions  'J',  =  a,  +  a2x  and  ft  =  ai  *  a2x  are  sufficiently 
differentiable,  independent,  complete,  and  values  for  the  coefficients  can  be 
found  to  satisfy  the  essential  boundary  conditions.  Below  the  Seredipity  method 
is  used  to  derive  the  interpolation  functions: 


The  boundary  conditions  are: 

+,(x=0)  =  1  t2(x=0)  =  0 

<f»i(x=h)  =  0  ^teh)  =  * 


Solving  for  coefficients: 
fi(0)  =  a,  =  1  *2(h)  =  a,  =  0 

«j»,(h)  =  1  +  a2h  =  0  ^(h)  =  a2h  =  1 


These  interpolation  functions  are  used  in  the  above  expressions  for  Mjj,  Kjj,  and 

Fj  to  derive  the  element  matrices.  In  the  derivation,  the  cross  section  of  the 

element  and  the  distributed  force  F(x,t),  are  allowed  to  vary  linearly  with  x.  The 
modulus  of  elasticity,  was  assumed  to  be  constant. 

[Kl  =  —  (A,  +  A2)  f  1  1  1 

2h  1-1  1  J 

_  T  (3m,  +  m2)  (m,  +  m2)  1 

12  L  (m,  +  m2)  (m,  +  3m2)  J 


Where  the  subscripts  indicate  values  at  the  left  and  right  end  of  the  bar  element. 
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The  governing  equation  for  the  beam  element  is 


82u  82  82u 

m  —  +  —  El  —  ♦  F(x,t)  -  0 

8t2  8x2  8x2 


Where  m  is  the  mass  per  unit  length,  F(x,t)  is  a  distributed  transverse  forcing 
functioa  Integrating  the  governing  equation  against  a  test  function  over  the 
domain  of  the  beam  element  gives 


82u  82  f  82u  1  1 

- fi  —  +  Ffxrt 


v  m— - -  El  — - 

8t2  8x2  8x2 


F(x,t)  d*  =  0 


Integrating: 

,h  , 


8v  8  82u  8  82u 

- El -  ♦vF(x,t)  dx*v —  El - 

8x  8x  8x2  8x  8x2 


The  last  term  of  the  above  expression  corresponds  to  the  shear  (natural  boundary 
condition)  at  either  end  of  the  element  and  will  be  denoted  as  Qi  (shear  at  the 
left  end)  and  Q3  (shear  at  the  right  end). 

Integrating  the  second  term  again: 


h 

_ 

x=h 

82u  82v  82u 

8v  82u 

vm  „  *  El  *vF(x,t) 

dx+vQ  j 

-El - - 

3t2  3x»  3x2 

i 

8x  8x2 

0 

■ 

x=0 

The  last  term  of  the  above  expression  corresponds  to  the  moment  at  either  end 
of  the  element  (natural  boundary  condition),  and  will  be  denoted  as  Q2  (moment 
at  the  left  end)  and  Q«  (moment  at  the  right  end). 
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The  displacement  is  again  interpolated  by  an  expression  of  the  form: 

2 

u  =  I  ui(t)  *j(x) 
j=1 

Substituting  the  above  expression  for  u,  and  v=tj(x)  results  in  the  matrix 
formulation 


[M]{u’  >  *  [KHuJ  =  (F(t)} 


Where  (’)  means  differentiation  with  respect  to  t  and 
.h 

M  lfl  .  iff  .  s4lf  1/  .  .  _ _  _ 

i>'J0  *»  dX* 


rn 

Mjj  =  J  m  fyfj  dx 


fh  d^j  d^j 

Kii  =  J  AE— 1  dx 


r" 

=  ’f'jFCx.t)  dx  ♦  Qj(t) 

J  o 


The  interpolation  functions  <J»j  (tor  i=l  to  4)  must  be  sufficiently  differentiable, 

independent  of  one  another,  complete,  and  must  satisfy  the  essential  boundary 
conditions.  The  expressions  tpai+a^+asx2*^*3  and  ^2=ai+a2x+a3x2+a4x3 
are  sufficiently  differentiable,  independent,  complete,  and  values  for  the 
coefficients  can  be  found  to  satisfy  the  essential  boundary  conditions.  Below 
the  Seredipity  method  is  used  to  derive  the  interpolation  functions. 


The  boundary  conditions  ( ’  denotes  differentiation  with  respect  to  x): 


♦iOpO)  =  1 
Mm*0)  =  0 
<h(x=h)  =  0 

«JV(x=h)  =  0 

Solving  for  the  coefficients: 

*,(0)  =  a,  =  I 

ti'(O)  =  a2  =  0 

fi(h)  =  I  ♦  a3h2  ♦  a4h3  =  0 

♦i’(h)  =  2a3h  ♦  3a4h2  =  0 


*2(jfO)»  0 
+2'(x=0)  =  -I 
t2(x=b)=  0 
t2’(x=h)  =  0 


♦z(0)  =  a,  =  0 

ta'CO)  =  a2  =-l 

f2  (h)  =  -h  ♦  a3h2  ♦  a4h3  =  0 

♦2'(h)  =  -1  ♦  2a3h  ♦  3a4h2  =  0 
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The  boundary  conditions  ( '  denotes  differentiation  with  respect  to  x): 


^3(x=0)  =  0 
♦stlFO)  =  0 
<j>3(x=h)  =  1 
*3*(x=h)  =  0 


^4  (x=0) 
M*=0) 

^4(x=h)  : 

M**) 


0 

0 

0 

-1 


Solving  for  the  coefficients: 
*3(0)  =  a,  =  0 
^3'(0)  =  a2  =  0 
hft)  =  a3h2  +  a4h3  =  1 
♦3‘(h)  =  2a3h  ♦  3a4h2  =  0 


♦4(0)  =  a,  =  0 
MO)  =  a2  =  0 
^4  (h)  =  a3h2  +  a4h3  =  0 
Mh)  =  2a3h  ♦  3a4h2  =  -1 


fh2 

h3  ■ 

a3 

[.  1 

rh2 

h3  1 

a3 

•  « 

0 

* 

,  =  A 

► 

(  * 

►  =  < 

bh 

3h2 . 

a4 

>  ■ 

l°J 

L2h 

3h2  J 

34 

■  « 

-1 

x3 

h2 


These  interpolation  functions  are  used  in  the  above  expressions  for  My ,  Kjj ,  and 

Fy  to  derive  the  beam  element  matrices.  The  cross  section  of  the  element  and 

the  transverse  loading  function  F(x,t)  are  allowed  to  vary  linearly,  but  the 
modulus  of  elasticity  is  held  constant. 


Page  *7 
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IK1  =  — 

h3 


IM1  =  — 
830 


6(1,  ♦  l2) 


(sym) 


-h(4l,  ♦  2I2)  -6(1,  ♦  12) 

h2(31,  +12)  h(41,  ♦  212) 

6(1,  ♦  I2) 


-h(2I,  ♦  4I2) 
h2(l,  ♦  l2) 
h(2l,  *  412) 
h2(I,  ♦  312) 


h(10m,t3m2)  -h2(151,  +  712)  h(9m,  ♦  9m2)  h2(7m,  ♦  6m2) 

h3(5m,  +3m2)  -h2(6m,  ♦  7m2)  -h3(m,  +  m2) 
(sym)  h(3m,  +  10m2)  h2(7m,  ♦  15m2) 

h3(3m,  ♦  5m2) 


h 

{F}  =  — 
60 


15(f,-3f2) 
-h(3f,  ♦  2f2) 
3(3f,  ♦  7f2) 
h(2f,  ♦  3f2) 


The  bar  and  beam  elements  are  now  superimposed  upon  one  another  to  form  the 
frame  element.  It  is  assumed  that  forces  and  displacements  in  the  axial 
direction  and  in  the  transverse  direction  are  independent  of  one  another.  The 
resulting  element  matrices  are  shown  below. 


h2(A,+A2) 


-h2(Ai+Ao) 


IK1  =  — 
2h3 


m  ,*l2)  -2W4I,+2I2) 

2h2(3l,+l2) 


-12(1, +l2)  -2h(2l,+4l2) 

2h(4l,+2l2)  2h2(l,+l2) 


h2 (A,+A2) 


(sym) 


12(1  ,+l2) 


2M2I,+4I2) 

2h2(li+lo) 
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ini= — 
840 


’70h(3mj+m2)  yOhtmj+mj) 

24h(10m1+3m2)  -2h^(15m  ,+7^)  54h(m  j+rr^)  2h^(7mj+6m2) 

h^(Smj*3m2)  -2h^(6m^+7m2)  -3h^(m,+fri2) 


Cf)  =  — 
60 


(sym) 

» 

Jon(2f,  ♦  f2) 
i5n(r,*3f2) 
-h2(3f,  ♦  2f2) 
10h(f,  ♦  2f2) 
3n(3f,  ♦  7f2) 
h2(2f,  ♦  3f2) 


70h(m^m2) 


24h(3m1-H0m2)  2h2(7m1+lSm2) 
h'H3m,+5m2) 


Prior  to  assembling  the  element  matrices  into  the  global  matrices,  the  element 
local  coordinates  must  be  converted  to  global  coordinates.  This  is  done  by 
premultiplying  the  stiffness  and  mass  matrices  with  the  following 
transformation  matrix.  The  element  force  matrix  is  premultiplied  by  the 
transpose  of  the  transformation  matrix.  The  angle  9  is  measured  from  the  from 
the  positive  x  direction  clockwise: 

cose  -sine  1 

cose 


(sym) 


cose  -sine 
cose 


When  the  element  matrices  are  assembled,  the  internal  element  forces,  P,,  Q,,  Q2, 
P2,  Q3,  and  Q«  are  canceled  out  by  the  internal  element  forces  of  adjoining 
elements.  There  may,  however,  be  externally  applied  loads  at  the  nodes,  if  this 
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is  the  case,  they  are  added  into  the  formulation  as  shown  in  the  above  expression 
for  the  force  matrix.  Note  however,  that  since  the  loads  are  applied  directly  to 
the  nodes,  that  the  coordinate  transformation  is  not  appropriate. 

DDlginq  Essential  Boundary  Conditions 
In  the  development  of  the  mass  and  stiffness  matrices,  the  shape  functions  were 
developed  in  order  to  account  for  essential  boundary  conditions  but  essential 
boundary  conditions  were  never  actually  applied.  As  a  result,  the  stiffness 
matrix  is  currently  singular  and  can  not  be  inverted  (ie,  the  problem  can  not  be 
solved  as  is).  One  consequence  of  this  is  that  this  configuration  can  not  be  used 
to  solve  for  displacements  of  structures  which  are  not  anchored  in  some  way  to 
an  immovable  object  (as  an  example  an  object  floating  in  space).  Application  of 
essential  boundary  conditions  constrain  the  structure  and  the  stiffness  matrix 
becomes  non-singular. 

There  are  two  approaches  to  apply  the  essential  boundary  conditions.  Since  the 
displacement  of  a  node  in  a  particular  degree  of  freedom  is  known,  the 
corresponding  equation  in  the  matrix  formulation  is  simply  changed  to  reflect 
the  known  value.  The  Guass  elimination  scheme  used  to  solve  the  simultaneous 
equations  will  insure  that  the  influence  of  the  displaced  node  is  properly 
reflected  thorough  out  the  structure.  This  is  the  method  used  in  the  program 
DynFEP.  It  has  the  advantage  that  all  constrained  nodes  need  not  all  move  at 
once  or  in  the  same  directions,  in  addition  rotations  of  individual  nodes  can  be 
investigated  with  this  approach. 

An  alternate  approach  is  described  in  the  referenced  text  by  Clough  &  Penzien. 
The  common  approach  used  in  earthquake  analysis,  is  to  drop  the  row  and  column 
corresponding  to  the  displaced/constrained  node  from  the  formulation,  reducing 
the  number  of  simultaneous  equations  to  be  solved.  Then  effects  of  base  motion 
are  added  into  the  formulation.  This  is  done  by  adopting  a  coordinate  system 
where  the  unknown  displacements  are  relative  to  the  movement  of  the  base  of 
the  structure.  Then  an  inertial  term  is  added  to  the  right  hand  side  of 
appropriate  equations.  As  an  example,  if  the  base  of  the  structure  experienced  a 
horizontal  displacement,  then  an  inertia  term  would  be  added  to  every  equation 
in  the  matrix  formulation  which  pertained  to  horizontal  displacements.  In 
matrix  formulation  an  acceleration  vector  accounting  for  horizontal  and  vertical 
movement  is  developed  and  premultiplied  by  the  mass  matrix  to  obtain  the 
inertia  term,  this  column  matrix  is  then  added  to  the  right  hand  side  of  the 
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equations. 

Rotations  are  normally  disregarded  in  this  approach.  First  because  earthquakes 
seldom  display  any  rotational  components  and  second  because  the  bookkeeping 
chore  is  very  burdensome.  The  inertia  effect  of  a  node  rotation  on  another  node 
is  proportional  to  the  lever  arm  between  the  two  nodes.  Thus,  for  each  node  that 
rotates  the  lever  arm  between  it  and  all  other  nodes  must  be  found  in  calculating 
the  inertia  term.  In  addition,  its  very  difficult  to  conceptualize  the  inertial 
effects  of  one  node  on  another  when  several  nodes  are  rotating. 

The  DynFEP.uncouple/solve  and  DynFEP.reduce  programs  presented  below  are 
formulated  in  the  above  manner. 
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Time  Approximations 

The  Finite  Element  Method  provides  a  method  of  converting  the  differentials  with 
respect  to  x  in  the  governing  equations  into  a  linear  algebra  problem  suitable  for 
solution  by  computer.  During  the  derivation  it  was  assumed  that  displacements 
with  respect  to  space  and  time  could  be  separated.  We  are  now  faced  with 
solving  the  resulting  matrix  differential  equation  in  time. 

lM]{u"}  *  (K]{u)  =  {F(t)> 

In  order  to  utilize  a  computer  based  solution,  the  above  differential  equation 
must  also  be  simplified  to  an  algebraic  form.  The  Newmark  method  of  direct 
integration  is  a  commonly  used  technique  to  accomplish  this.  The  Newmark 
method  is  described  in  referenced  texts  by  Reddy.  Clough  &  Penzien,  and  Bathe  & 
Wilsoa  It  is  based  on  the  following  assumptions: 

(U*>t+At  =  (u‘}t  ♦  1(1  -  8){u*'}t  ♦  8{u%At]At  (I) 

{u>t+At  =  {u)t  ♦  (u'}tAt  ♦  [i}/2  -  <x){u")t  ♦  <x{u"}t+At)At2  (2) 

Where  <x  and  8  are  parameters  that  can  control  the  integration  accuracy  and 
stability.  When  8=,/2  and  <x=Vs  the  above  expressions  correspond  to  a  linear 
acceleration  assumption.  When  S=,/2  and  <x=,/4  above  expressions  correspond  to 
a  constant-average-acceleration  assumptioa 

Working  with  equation  (2),  acceleration  for  a  new  time  increment  can  be 
expressed  in  terms  of  current  displacement  and  values  from  the  last  time 
increment. 

oc(u,,}t4  tAt  =  (u)t+At  -  (u}t  -  {u’jAt  -  (V2  -  <x){u")t  At 


Cu)t)  -  a2{u*}t  -  a3{u")t 


•  i  - 
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Substituting  equation  (3)  into  the  discretized  equations  of  motion: 

[Ml(a1({u)t+At  -  «ult)  -  a2{u*)t  -  a3{u")t)  ♦  lKKu)t+Al  =  {F)t+At 

(aiini  +  lKl){u)t+;t  =  (F)t+At  ♦  lM]({u}t  ♦  a2{u')t  +  a3{u''}t)  (4) 

Using  the  above  equation  the  procedure  for  direct  integration  is  as  follows: 

1)  Knowing  displacement,  velocity,  and  acceleration  from  the  last  time 
step  (or  from  initial  conditions),  find  displacements  for  next  time  step  using 
equation  (4)  above. 

2)  Using  equation  (3)  find  current  acceleration. 

3)  Using  equation  (1)  find  current  velocity. 

4)  Proceed  to  next  time  step. 

The  Newmark  method  is  unconditionally  stable  for  <x=V2  and  S2’/*  and  is 
normally  stable  for  <x=V2  and  S=V6.  In  order  to  also  insure  accuracy  of  the 
method.  At  should  not  exceed: 
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Modal  Analysis 

The  Finite  Element  Method  and  the  Newmark  method,  are  used  above  to  convert 
the  differentials  which  govern  movement  of  plane  frame  structures  to  a  set  of 
simultaneous  algebraic  equations.  These  equations  are  then  solved  repeatedly  in 
small  time  steps  to  obtain  the  displacement  response  of  the  structure  over  time. 
Given  this  method  of  solution  it  should  be  obvious  that  any  means  to  further 
simplify  the  solution  process  will  be  valuable. 

The  texts  by  Clough  &  Penzien,  and  Bathe  &  Wilson  present  a  widely  used  method 
to  uncouple  the  simultaneous  equations  so  that  they  may  be  solved  independently 
of  one  another.  The  method  involves  expressing  the  equations  of  motion  as  an 
eigenvalue  problem,  solving  the  eigenvalue  problem,  and  then  re-expressing  the 
equations  of  motion  in  a  coordinate  system  which  has  been  generalized  by  the 
eigen  vectors. 


If  the  structure  in  question  is  in  free  vibration  then  the  forces  on  the  right  hand 
side  of  the  equations  of  motion  are  equal  to  zero,  [MHu”MK]{u}={0}.  The 
solution  for  each  degree  of  freedom  is  then  {u}={<p}sin(a)t).  Substituting  this 
solution  into  the  equations  of  motion 

-o)2  [M]{<p}sin(a>t)  ♦  lK]{<p}sin(<i)t)  =  {0} 

([K]  -a>2  lMl){<p}sin(<i>t)  =  {0} 

Since  sin(a)t)  is  not  equal  to  zero  for  all  t, 

«KJ  -o>2lM]){<p}  =  {0} 


A  non-trivial  solution  to  this  system  of  simultaneous  equations  exists  only 
when  |  |K]  -a>2[Ml|=0.  When  this  determinate  is  expanded,  it  results  in  an 
algebraic  equation  of  the  n&  degree  (where  the  dimensions  of  (K1  and  iMl  are 
n-by-n).  The  n  roots  to  this  equation,  o>j  (where  i=l,  2,  .  .  .,  n),  represent  the 

frequencies  of  the  n  modes  of  vibration  that  are  possible  in  the  system.  The 
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associated  eigen  vectors,  describe  the  relative  displacements  of  the 

structure  nodes  in  the  response  mode.  The  total  response  is  given  by  the  sum 
of  the  mode  responses  each  multiplied  by  a  currently  unknown  amplitude. 

The  eigen  vectors  are  [K]  and  [M]-orthogonal.  Thus  {<Pj}T[Ml{<pj}=[Mn],  where 
mnl  is  a  diagonal  matrix,  and  {<Pj}T[KH<Pj}=[Knl,  where  [Knl  is  a  diagonal  matrix. 
In  addition  {<Pj}T[M]{<pj}=l01,  and  {<Pj}T[K]{<Pj)=l01  where  i*j. 

The  advantage  of  the  modal  analysis  is  seen  when  a  generalized  coordinate 
system  is  defined  as  {u}=l*m)-  Where  [$]  is  a  matrix  made  up  of  the  individual 
eigen  vectors.  Premultiplying  the  original  equations  of  motion  by  [*1T  and 
substituting  the  generalized  coordinate  system  into  the  equation  of  motion 
results  in: 

♦  MT[KimU)  =  MT{F(0) 

[MnlU ")  *  (KnM?)  =  Mt(F(0) 

Stated  in  terms  of  the  generalized  coordinate  system,  the  equations  of  motion, 
are  uncoupled.  Since  [Mnl  and  [Kn]  are  diagonal  matrices  each  equation  in  the 

above  system  of  equations  is  independent  of  the  others. 

Solution  of  the  Eigenvalue  Problem 

Clough  &  Penzien  describe  a  matrix  iteration  method  originally  developed  by 
Stodola  to  solve  the  eigenvalue  problem.  The  eigenvalue  problem  is  restated  as 
follows: 

KM?)  =  <*)2[nK<p> 

Rearranging: 

WlMM?)  =—,(?> 

0>2 


The  Stodola  method  consists  of  using  a  guessed  trial  mode  shape,  C^Ptrial  on 
the  left-hand  side  of  the  above  equation  to  calculate  a  new  guess  on  the 
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right-hand  side.  The  square  or  the  frequency  is  obtained  by  dividing  any 
component  of  the  new  guess  by  the  same  component  of  the  original  guess.  The 
new  guess  will  always  be  better  then  the  old  guess,  and  the  process  will 
converge  to  the  lowest  mode  or  frequency. 

Using  the  orthogonal  properties  of  the  eigenvectors  it  is  possible  to  eliminate 
the  components  of  any  particular  mode  from  the  total  response  of  the  structure. 
By  eliminating  the  first  mode  components  from  the  total  response  it  is  possible 
to  use  the  above  method  to  find  the  second  mode  response  (since  it  v'ould  now  be 
the  lowest).  Extending  this  approach,  succeeding  modes  can  also  be  found. 

Expressing  a  trial  mode  shape  in  terms  of  its  modal  components  and  then 
premultiplying  both  sides  by 

n 

Wtrial  >  =  Z  WjJAj  =  {<p,}A,  ♦  {<p2)A2  ♦  I«P3)A3  ♦  .  ♦  {^Ap 

i=l 

{<MT[nl{<ptr  ja,  }=  {ty}TlM]{<p,}A,*{<p,)%il{<p2}A2*  •  •  ♦  (f1}T[nHfn)An 


{<Pt}TlMl{<Ptria,  }=  1<P|)TIM]{<P,}A, 

Solving  for  A,: 

A  _  W^tiHftrial  > 

1 "  {'P|}TIM1{<P,J 

Subtracting  the  first  mode  shape  from  the  original  trial  mode  shape  results  in  a 
new  trial  with  no  first  mode  components,  {<p  trials  ^ 

{<?,){?, )T[nl 

Wtrialm 1  =  Atrial  >  ‘  WAi =  Atrial 1 "  ^ial  > 


This  can  also  be  expressed  as  C^Ptriald,  )  =  [SiH'Ptrial  wher* 


IS,]  =  II]  - 


(jWiTjni 
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The  (S11  matrix  is  referred  to  as  the  first  mode  sweeping  matrix.  It  has  the 
property  that  when  multiplied  by  any  trial  vector  it  removes  the  first-mode 
component.  Sweeping  matrices  which  remove  more  then  one  mode  shape  can  be 
constructed  in  a  similar  matter.  As  an  example,  a  sweeping  matrix  which  will 
remove  the  first,  second,  and  third  mode  shape  components  from  a  trial  vector 
would  be  constructed  as  follows: 

. .  WiM,MTtnl  C«p2H<p2lT[Ml 


W'lMlW,}  {<P2>t[mH<p2}  {<pj>TtM|C<p3> 

The  resulting  Stodola  matrix  iteration  model  to  find  the  fourth  mode  shape  and 
eigenvalue  becomes: 

IKI-'IMHSjIW)  =—„<?} 


The  method  of  solution  suggested  by  the  above  methods  consists  of  the 
following: 

1)  Find  lowest  mode  shape  and  eigenvalue  using  the  matrix  iteration. 

2)  Using  the  newly  calculated  first  mode  shape  eliminate  the  first  mode 
components. 

3)  Repeat  the  procedure  for  the  next  mode  shape  and  eigenvalue. 

Each  successive  mode  shape  is  based  on  eliminating  the  previous  mode’s 
components.  According  to  Clough  &  Penzien,  numerical  roundoff  errors  which 
allow  any  previous  mode  components  to  remain  in  the  sweeping  matrix  are 
accumulative.  Thus  In  order  for  the  sweeping  matrix  to  perform  effectively  for 
higher  modes  it  Is  necessary  to  retain  a  great  deal  of  precision  in  calculating  the 
lower  modes. 


The  eigenvalues  and  eigen  vectors  are  now  used  to  form  the  uncoupled  equations 
of  motion.  The  Newmark  method  is  applied  to  the  resulting  independent 
equations.  The  independent  equations  are  solved  in  the  terms  of  the  generalized 
coordinates,  {£}.  while  stepping  through  time,  in  each  time  step  the  real 
displacement  vectors  are  found  from  the  relation  {u>  =  [#]{£}• 
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Reduction  of  the  Equations  of  Motion 

When  the  above  methods  are  applied  to  real  structures,  very  large  matrices  and 
correspondingly  large  computer  capacity  are  required  to  solve  the  resulting 
equations.  It  is  thus  desirable  to  reduce  the  number  of  equations  which  must  be 
solved.  In  the  study  of  structures  it  has  been  determined  that  only  the  first  few 
response  modes  contribute  significantly  to  the  overall  response  a  structure.  It 
is  therefore  reasonable  to  ignore  the  higher  modes  of  response  if  it  will  reduce 
the  number  of  equations  to  be  solved. 

Robert  J.  Guyan  described  such  a  method  of  reducing  the  number  of  equations  to 
be  solved  in  a  paper  to  the  AIAA  Journal.  The  method  consists  of  a  static 
reduction  of  the  structure.  Working  with  the  static  description  of  the  structure, 
the  matrices  are  partitioned  by  the  nodes  which  will  be  retained  in  the  solution 
(referred  to  as  the  primary  nodes),  and  the  nodes  which  will  be  eliminated  from 
the  formulation  (referred  to  as  the  secondary  nodes).  It  is  assumed  that  no 
external  loads  will  be  applied  to  the  secondary  nodes. 

T  Kppl  iKpsl  1  f  fup)  1  _  |  1 

L  [Kspl  I^ssl  ■*  |  <%>  J  1  J 

This  results  in  the  following  two  matrix  equations: 

[KppHUp}  +  iKpgHUg)  =  (fp) 

[KgpHUp}  ♦  fKggHUg)  =  {fg}  =  {0} 

Multiplying  the  second  equation  by  lKpSHKssr1 

lKpg][Kgg]  '  [Kgp]{Up)  +  [Kpg][Kgg]  '  [Kggl{Ug)  =  (0) 

[KpgHKgg]  '  IXgpIlUp!  +  CKpgHUg)  =  (0) 


(5) 

(6) 


(7) 
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Subtracting  equation  (7)  from  equation  (5)  gives: 

[K*J{up}  =  {fp} 

Where  IK*1  is  the  reduced  stiffness  matrix,  found  by  the  following  expression: 
IK-1  =  [Kppl  -  [KpS][Kssr'  lKspl 

In  addition,  from  equation  (7)  a  transformation  matrix  can  be  obtained  to 
convert  between  the  primary  and  secondary  displacement  values. 

[Kps^s^  '  -KpslKssl  '  ^sp^up^ 

(us)  =  -Kssr'  IKgpHUp}  =  -ITMUp)  (8) 

Rearranging: 


{u}  = 


III 

(Up) 

-IT] 

►  * 

{up} 

The  kinetic  energy  of  the  structure  can  be  expressed  as: 
T  =  V2{u'}t[MKu'} 

Substituting  the  above  expression  for  Cu}  results  ire 


T  =  VS 


111 

r  CMppi  iMpsi  l 

111 

-IT] 

L  l^spl  {Nsgl  ■* 

-1T1 

Cu  p) 
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Thus  it  can  be  seen  that  the  reduced  [Ml  can  be  expressed  as: 


Expanding  the  above  expression  and  substituting  the  expression  for  [Tl  we  obtain 
the  simplified  expression  for  [M*l : 

|M*1=  iMpp]-  [MpgllKgg]  '  iKgp]  -  [  KpgllKggl  ^  ([Mgpl~  [Mgg][Kggl  ^  [Kgpl)  (9) 

The  above  method  has  reduced  the  mass  and  stiffness  matrix  of  the  structure  and 
therefore  the  number  of  equations  which  must  be  solved.  However,  the 
transformation  matrix  [T]  used  to  find  the  displacement  of  secondary  nodes  has 
ignored  any  inertial  effects.  The  exact  expression  for  [Tl  is  found  by  expressing 
the  eigenvalue  problem  in  the  partitioned  form: 

[Kppl  [Kpsl  1  1 9P)  _  ^  [Mppl  [MpSl  j  {<pp}  _  ^ 

.  [Kgpl  [Kggl  J  (fg)  .  [Mgpl  [Mggl  J  (^g) 

Working  with  the  second  partitioned  matrix  equation  the  exact  transformation 
matrix  for  the  eigen  vectors  is  obtained: 

lKspH<pp>  *  IKSSJ(<PS}  -  oJmspH<Pp}  -  <i>2[Mss]{4>s}  =  {0} 

(<0*IMSS1  -  IKSS1)(<PS1  =  (d)2lnspl  -  IKSP])WP) 

IT]  =  (u2mss]  -  IKSS1)-'  (d)2[nSpl  -  IKspl) 

Note  that  if  the  inertial  terms  in  the  above  expression  are  neglected,  the  same 
transformation  matrix  developed  earlier,  based  on  a  static  derivation,  is 
obtained.  The  above  expression,  however,  involves  an  eigenvalue  based  on  the 
complete  set  of  equations  and  requires  that  an  inversion  of  the  (o>2[Mssl  -  [Kssl) 

term  be  found  for  each  eigenvalue. 
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Charles  Miller  describes  a  transformation  matrix  which  is  more  accurate  then 
the  one  derived  previously  and  more  convenient  then  the  exact  formulation  in  a 
paper  to  the  Journal  of  the  Structural  Division,  Proceedings  of  the  ASCE. 

Mr.  Miller  notes  that  0)2IMSS1  and  ti>2lMspI  are  normally  small  when  compared  to 
[Kssl  and  [Kgpl.  With  this  in  mind  he  expands  the  (a)2[MSg]-[KgSl)-1  term  of  the 
exact  formulation  about  iKgg]'1  dropping  the  terms  in  comparison  to  a>2 
terms.  This  results  in  a  revised  [T]: 

in  =  iKssr'  iKSpi  *  w2(-iKssr'  ttisp]  *  iKssr'  imssmkssi-1  [kspd 

Expressing  the  first  equation  of  the  partitioned  eigenvalue  problem  in  terms  of 

Wp) 

[KpplWp}  -  [KpgHTMp}  -  a)2lMppH<pp}  -  u)2[MpS][TK<Pp}  =  {0} 

Expanding  the  o)2lMpS][T]{<Pp}  term  again  dropping  the  terms 

^IMPSHTKV  =  "^2I^psHKggl  ^ 

Expanding  the  [KpSHTl{<pp)  term 

[KpsITMfp)  -  ([KpsI[Kssl  ^  ^sp1 

~&>2lKpgl(“[Kgg|  ^  [HgpMKgg]  ^  ^SS^SS^  ^  ^gplM^Pp) 

Substituting  these  expanded  expressions  into  the  eigenvalue  problem: 

(lKpp] +  [KpS][Kssl  '  iKgpl  "^[KpgK'lKgg]  ^  tMspMKggj  ^  iMgsHKgs^  '  ^spl 

-  “2i»pPi  *  <o2[tipS][Kssr'  uispm<Pp)  =  (oj 
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Rearranging  gives: 

([Kppl  *  iKps^ss^  ^  [KgplH^pJ  - 

^(iNpp^l^ps^ss^  ^pHSs^ssl  '  (iMgpl^lMgs^^ss^  1  iKspDC^Pp) 
iK**H<pp>  =  ^rH^pj 

The  revised  transformation  matrix  results  in  the  same  expressions  for  the 
reduced  mass  and  stiffness  matrices.  The  revised  transformation  matrix 
requires  only  one  inverse  of  a  partition  of  the  stiffness  matrix,  not  a  new 
inverse  for  each  eigenvalue.  In  addition,  the  revised  transformation  matrix 
should  provide  more  accurate  displacements  since  it  partly  accounts  for  inertia 
terms.  The  more  accurate  displacements  provide  a  more  accurate  basis  for 
approximating  internal  forces. 

The  process  suggested  by  the  above  formulation  proceeds  as  follows: 

1)  Partition  mass  and  stiffness  matrices  and  find  reduced  matrices  using 
equations  (8)  and  (9). 

2)  Solve  the  eigenvalue  problem  for  reduced  mass  and  stiffness  matrices. 

3)  For  each  eigenvalue  find  the  transformation  matrix. 

4)  Use  each  transformation  matrix  to  obtain  full  eigen  vector. 

Once  the  full  eigen  vectors  are  found  the  solution  proceeds  the  same  as  under 
Modal  Analysis.  It  should  be  noted,  however  that  the  dimensions  of  various 
matrices  have  been  changed. 

Where  n  equals  the  number  of  unknowns  in  the  structure,  and  m  equals  the  number 
of  modes  retained  in  the  solution,  the  dimensions  of  the  eigenvalue  matrix  is 
m*m  and  the  dimension  of  the  complete  eigen  vector  matrix  (*]  is  n*m.  When  the 
generalized  mass  matrix  is  found  from  the  relation  1*]TIM][*1,  its  dimensions 
are  mm  Due  to  the  orthogonality  of  the  eigen  vector  matrix,  the  generalized 
mass  matrix  is  still  diagonal  and  there  remain  only  m  independent  equations  to 
be  solved.  Converting  the  generalized  solutions  to  real  coordinates  using  the 
relation  {u}  =  1*1(0  results  in  full  size  displacement  matrix  ((*1  is  dimensioned 
n*m  and  (0  is  dimensioned  m*i). 
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The  Dynamic  Finite  Element  Program 

The  Dynamic  Finite  Element  Program  is  a  system  of  programs  developed  to  apply 
the  methods  presented  above.  The  programs  are  written  in  Micro-Soft  Basic  for 
the  Apple  Macintosh,  version  2  (Micro-Soft  Inc.  is  currently  developing  versions 
of  this  advanced  version  of  Basic  for  IBM  compatible  machines).  Flow  diagrams 
and  listings  of  the  programs  are  provided  in  Appendix  A.  A  description  of  each 
program  and  its  operation  follows. 


The  DynFEP.menu  program  serves  to  connect  the  system  of  programs  together.  It 
provides  a  menu  from  which  the  user  can  choose  to  create  data  files  which 
describe  structures,  or  to  solve  problems  which  have  been  defined  earlier. 
Problems  may  be  solved  in  any  of  three  ways,  using  dynamic  reduction,  using 
modal  analysis,  or  direct  numerical  integration  of  the  equations  of  motion  using 
the  Newmark  method. 

The  DynFEP.menu  program  maintains  control  of  the  flow  of  execution  by  passing 
five  variables  to  each  program  in  the  system.  The  five  variables  are:  the  number 
of  global  nodes  labeled  as  GN;  the  number  of  elements  labeled  as  NE;  the  number 
of  unknowns  remaining  after  application  of  the  essential  boundary  conditions 
labeled  as  n  (if  essential  boundary  conditions  are  applied  by  reducing  the  number 
of  equations);  the  number  of  modes  to  be  retained  in  the  solution  labeled  m  (if 
the  structure  is  to  be  reduced);  and  a  string  variable  describing  the  chosen 
solution  method  labeled  as  Path$. 


This  program  creates  the  data  files  which  describe  the  structure,  and  forcing 
functions  and  displacements  applied  against  it.  Data  describing  the  structure  is 
entered  using  Basic  DATA  statements.  A  separate  Basic  program  listing 
containing  the  only  the  desired  DATA  statements  is  prepared  and  saved  under 
ASCII  (text  only)  format.  The  program  assumes  that  such  a  program  has  been 
prepared  and  merges  with  it.  When  the  resulting  new  program  is  executed,  it  will 
read  the  prepared  data  and  create  the  required  structure  data  files. 
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The  DATA  statements  must  be  formatted  to  match  the  READ  statements  in  the 
>  DynFEP.create  program.  This  is  normally  accomplished  by  copying  a  previously 

created  set  of  DATA  statements,  and  modifying  them  to  fit  the  new  problem 
using  Basic’s  editing  capabilities. 

DynFEP.mass/stif  f  ness 

This  program  reads  previously  created  structure  data  files  and  assembles  the 
global  stiffness  matrix,  the  global  mass  matrix,  and  the  global  static  forces 
matrix.  The  program  steps  through  the  structure  elements,  and  constructs  the 
element  matrices  using  the  relations  presented  above.  The  orientation  of  the 
element  is  checked  and  the  matrices  are  transformed  if  necessary.  Then  the 
element  matrices  are  assembled  into  the  global  matrices  in  accordance  with 
their  end  node  points. 

The  program  will  be  executed  as  determined  necessary  by  the  DynFEP.menu 
program.  The  above  matrices  will  be  assembled  once  for  any  particular 
structure,  it  is  not  necessary  to  reassemble  the  global  matrices  for  different 
applied  dynamic  forces  or  specified  displacements.  When  the  program  completes 
assemblage  of  the  global  matrices  it  will  check  a  specified  solution-pathway 
which  was  set  by  the  DynFEP.menu  program.  There  are  two  possible  paths, 
directly  to  DynFEP  if  non-modal  analysis  is  to  be  done,  or  to  DynFEP.essential  BC 
if  modal  analysis  is  to  be  done.  The  program  will  chain  to  the  appropriate 
program. 

DynEEP.essential  BC 

This  program  reads  the  structure  node  information  file,  determines  where 
essential  boundary  conditions  are  to  be  applied,  and  then  applies  the  conditions 
by  eliminating  the  appropriate  rows  and  columns  of  the  global  mass,  stiffness, 
and  static  force  matrices.  The  program  also  creates  a  boundary  condition  index 
which  will  be  used  by  succeeding  programs  to  reduce  the  global  dynamic  force 
matrix,  and  to  add  the  inertial  effects  of  moving  nodes  into  the  global 
formulation  (ie,  to  finish  applying  the  essential  boundary  conditions). 

The  program  will  be  executed  if  the  chosen  solution  method  is  modal  analysis  or 
dynamic  reduction  of  the  structure.  When  the  program  completes  its  execution, 
it  will  check  the  specified  solution-pathway  and  chain  to  the  appropriate 
program.  There  are  two  pathways  possible,  the  program  will  chain  to 
DynFEP.eigen  solver  if  modal  analysis  is  the  chosen  solution  method  or  it  will 
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chain  to  DynFEPxeduce  if  dynamic  reduction  is  the  chosen  solution  method. 
DynFEP.reduce 

This  program  reads  the  reduce  index  created  by  DynFEP.create  (from  user  input) 
and  reduces  the  number  of  equations  using  the  methods  presented.  In  addition 
the  program  prepares  two  matrices  (stored  in  temporary  disk  files)  which  are 
used  by  DynFEP.uncouple/solve  to  transform  the  primary  eigen  vectors  (eigen 
vector  of  the  reduced  structure)  into  a  eigen  vector  describing  the  full  structure. 

The  program  is  executed  along  the  dynamic  reduction  solution  pathway.  It 
executes  after  DynFEP.essential  BC.  When  it  completes  execution  it  chains  to  the 
DynFEP.eigen  solver  program. 

DynFEP.eigen  solver 

This  program  loads  prepared  mass  and  stiffness  matrices  and  solves  the 
corresponding  eigen  value  problem  using  the  Stodola  method  and  a  sweeping 
matrix  as  presented  above.  The  result  of  the  program  are  competed  matrices  of 
the  eigenvalues  and  eigen  vectors. 

The  program  assumes  that  the  global  matrices  have  had  the  rows  and  columns  of 
constrained  degrees  of  freedom  (ie,  specified  static  and/or  dynamic 
displacement)  removed.  The  program  uses  the  five  variables  passed  to  it  by  the 
DynFEP.menu  program  to  determine  if  the  global  matrices  have  been  reduced.  If 
the  matrices  have  been  reduced,  the  program  finds  the  transformation  matrix  for 
each  mode  frequency  (using  information  prepared  by  DynFEP.reduce)  and  uses  it 
to  transform  the  partial  eigen  vectors  into  a  full  eigen  vectors. 

The  DynFEP.eigen  solver  will  execute  if  the  dynamic  reduction  or  modal  analysis 
method  of  solution  is  chosen.  The  program  chains  to  DynFEP.uncouple/solve 
upon  completion. 

DynFEP.uncouple/solve 

This  program  pulls  together  the  work  of  previous  programs  and  solves  the 
problem  to  its  completion.  Its  execution  results  in  a  displacement  -vs-  time 
history,  for  each  node,  stored  on  disk.  The  displacements  recorded  are  relative 
to  the  movement  of  the  base.  Information  provided  by  the  user  regarding 
specified  note  displacement  is  assumed  to  be  accelerations  of  uniform  base 
movement. 
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The  DynFEP.uncouple/solve  program  loads  the  static  force  matrix,  the  boundary 
conditions  index,  the  mode  shape  matrix,  the  eigenvalues,  and  the  mass  matrix. 
If  the  structure  has  not  been  reduced  the  program  also  loads  the  initial 
conditions  and  converts  them  to  the  generalized  coordinates  to  be  used  to  start 
the  Newmark  direct  integration  scheme  (If  the  structure  is  reduced  the  mode 
shape  matrix  is  not  square  and  can  not  be  inverted  to  find  the  initial  conditions 
generalized  form.  Therefore,  if  the  structure  is  reduced,  all  initial  conditions 
must  equal  zero).  The  program  then  finds  the  generalized  mass  matrix  and  starts 
the  numerical  integration  scheme. 

The  first  operation  for  each  time  step  is  to  find  the  new  dynamic  force  matrix 
including  the  inertial  effects  of  base  motion.  The  dynamic  forces  for  applied  to 
nodes  or  members  are  first  found  for  every  node  in  the  structure;  the  dimension 
of  the  matrix  in  this  form  is  3(GN)*I  (where  3(GN)  means  three  times  the  number 
of  global  nodes).  Then  the  essential  boundary  conditions  are  applied,  resulting  in 
a  n*1  matrix  (where  n  is  equal  to  the  number  of  unknown  node  displacements). 
The  inertia  forces  are  obtained  by  multiplying  a  base  acceleration  matrix  by  the 
mass  matrix  (after  essential  boundary  conditions,  therefore  matrix  is  n*l).  The 
dynamic,  inertial,  and  static  force  matrices  are  then  summed.  Finally  the  force 
matrix  is  transformed  to  its  generalized  form  by  premultiplying  it  with  the 
transpose  of  the  mode  shape  matrix.  The  dimension  of  the  generalized  force 
matrix  is  either  n*1,  or  m*l  if  the  structure  has  been  reduced  (where  m  is  equal 
to  the  number  of  retained  modes). 

The  equations  are  now  in  their  uncoupled  form.  The  Newmark  method  is  applied 
to  solve  for  the  generalized  displacements,  velocities,  and  accelerations  for  the 
current  time  step.  Once  the  generalized  displacements  are  calculated  the  real 
displacements  are  found  by  premultiplying  the  generalized  displacements,  and 
their  derivatives,  by  the  mode  shape  matrix.  The  real  displacements  are  then 
stored  and  the  program  goes  to  the  next  time  step.  The  program  proceeds  for  a 
specified  number  of  time  steps. 

The  DynFEP.uncouple/solve  program  is  the  ending  program  for  either  the  dynamic 
reduction  solution  method  or  modal  analysis.  Upon  completion  the  program 
chains  to  the  DynFEP.menu  program. 
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DynFEP 

The  DynFEP  program  loads  the  full  mass  and  stiffness  matrices,  and  the  initial 
conditions.  Using  the  Newmark  method  and  Guass  elimination  it  solves  the 
problem  in  its  complete  form.  The  result  of  the  program  is  a  displacement  -vs- 
time  history  of  every  node  in  the  system.  The  displacements  recorded  are 
absolute  with  regards  to  the  coordinate  system.  Information  provided  by  the 
user  regarding  movement  of  nodes  is  assumed  to  absolute  displacement  also. 

The  first  operation  for  each  time  step  is  to  find  the  new  dynamic  force  matrix. 
The  program  reviews  the  node  and  element  loading  and  displacement  information 
stored  in  the  structure  data  files  (displacement  information  is  assumed  to  be 
absolute  displacement).  If  appropriate  time  history  files  for  non-harmonic 
forces  will  also  be  accessed.  The  new  dynamic  force  matrix  is  constructed  and 
added  to  the  static  force  matrix.  During  this  process  the  program  also 
constructs  a  boundary  condition  index. 

The  Newmark  method  is  applied,  then  the  boundary  condition  index  is  used  to 
apply  the  essential  boundary  conditions.  The  essential  boundary  conditions  are 
applied  by  modifying  the  equations  which  express  the  known  displacement.  The 
program  will  again  access  node  information  stored  on  disk  to  find  the  specified 
displacement,  accessing  time  history  files  where  appropriate  for  non-harmonic 
displacement  of  nodes. 

The  equations  are  now  solved  using  a  Guass  elimination  technique,  time  is 
incremented  and  the  process  is  repeated.  The  program  continues  for  a  specified 
number  of  time  steps. 

The  DynFEP  program  executes  ority  when  this  method  of  solution  has  been  chosen. 
When  the  program  completes  execution,  it  returns  control  to  the  DynFEP.menu 
program. 

Daia  Files 

Information  describing  the  structure  are  contained  in  five  primary  files,  the 
information  file,  the  node  file,  and  the  element  file,  initial  conditions,  and 
reductions  (when  dynamic  reduction  is  to  be  used).  These  files  are  created  by 
DynFEP.create  from  information  provided  by  the  user.  A  summary  of  these  data 
files  and  their  structure  is  presented  below. 
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Information  Data  File: 

Name:  structure  Name> 

Type:  permanent,  sequential  text 


Info: 

Descriotion 

Field 

Length 

Field 

Ngm? 

Remarks 

number  of  global  nodes 

n/a 

number  of  elements 

n/a 

number  unknown  displacements  n/a 

number  of  retained  modes 

n/a 

Node  Data  File: 

Name:  structure  Name>. Nodes 

Type:  permanent,  random  access 

Info: 

Field 

Field 

DescriDtion 

Lenqth 

Name 

Remarks 

Flag  describing  boundary  cond. 

12 

Fig« 

string  variable 

X  global  coordinate 

4 

N$0) 

single  precision 

Y  global  coordinate 

4 

N$(2) 

single  precision 

Horz.  Static  Load  or  Displ. 

4 

N$(3) 

single  precision 

Dyn.  Load  or  Displ.  Amp. 

4 

N$(4) 

single  precision 

dynamic  frequency 

4 

N«5) 

single  precision 

dynamic  phase  angle 

4 

N$(6) 

single  precision 

time  history  file  name 

8 

N$(7) 

string  variable 

Vert.  Static  Load  or  Displ. 

4 

N«8) 

single  precision 

Dyn.  Load  or  Displ.  Amp. 

4 

N$(9) 

single  precision 

dynamic  frequency 

4 

N$(tO) 

single  precision 

dynamic  phase  angle 

4 

N$(11) 

single  precision 

time  history  file  name 

8 

N$(t2) 

string  variable 

Rot.  Static  Load  or  Displ. 

4 

N«13) 

single  precision 

Dyn.  Load  or  Displ.  Amp. 

4 

N$(14) 

single  precision 

dynamic  frequency 

4 

N«1S) 

single  precision 

dynamic  phase  angle 

4 

N$(16) 

single  precision 

time  history  file  name 

8 

N$(17) 

string  variable 
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Meaning  of  flag  variable: 


ttt 


- Horiz.  load  or  displ.  (1=load,  0=displ.) 

- Horiz.  static  load  or  displ.  (1=yes,  0=no) 

- Horiz.  dynamic  load  or  displ.  (!=yes,  0=no) 

- Is  dynamic  load  harmonic?  (I=yes,  0=no) 

- Vert,  load  or  displ.  (l=load,  0=displ.) 

- Vert,  static  load  or  displ.  (1=yes,  0=no) 

- Vert,  dynamic  load  or  displ.  (1=yes,  0=no) 

- Is  dynamic  load  harmonic?  (1=yes,  0=no) 

- Rotational  load  or  displ.  (1=load,  0=displ.) 

- Rotational  static  load  or  displ.  (I=yes,  0=no) 

—  Rotational  dynamic  load  or  displ.  (I=yes,  0=no) 
r  Is  dynamic  load  harmonic?  (l=yes,  0=no) 

_ 12  character  flag  (string  variable). 


The  above  file  structure  allows  the  user  to  specify  both  a  static  and  a  dynamic 
load  or  displacement  at  any  node  (the  modal  methods  of  solution  do  not  support 
independent  displacement  of  nodes).  The  DynFEP  programs  will  interpret  the 
stored  information  to  be  either  a  specified  load  or  a  specified  displacement 
depending  on  the  above  twelve  character  flag.  The  flag  also  tells  the  program  to 
whether  or  not  to  look  for  static  or  dynamic  loads  and  whether  the  dynamic 
loads  are  harmonic  or  non-harmonic.  The  inclusion  of  a  phase  angle  allows 
nodes  to  be  loaded  or  displaced  out  of  phase  of  one  another  for  harmonic 
displacement  or  loading. 

Element  Data  File: 

Name:  structure  Name>.EIements 

Type:  random  access 


Info: 

Description 

Field 

Length 

Field 

Remarks 

Flag  describing  element  loading 

6 

Flg2$ 

string  variable 

Global  node  number  of  left  side 

2 

Lt$ 

integer 

Global  node  number  of  right  side 

2 

Rt$ 

integer 

Page  *29 


CE-685 
Term  Project 


Larry  Goshorn 
August  1985 


Descriotion 

Field  Field 
Length  Name 

Remarks 

Left  side  moment  of  inertia,  It 

4 

E«1) 

single  precision 

cross  sectional  area,  A, 

4 

E«2) 

single  precision 

mass  per  unit  length,  m5 

4 

E«3) 

single  precision 

Right  side  moment  of  inertia,  I2 

4 

E$(4) 

single  precision 

cross  sectional  area.  A2 

4 

E$(5) 

single  precision 

mass  per  unit  length,  m2 

4 

E$(6) 

single  precision 

Modulus  of  elasticity,  E 

4 

E$(7) 

single  precision 

Transverse  static  load  left  side 

4 

E$(8) 

single  precision 

static  load  right  side 

4 

E$(9) 

single  precision 

dynamic  amplitude 

'  4 

E$(I0) 

single  precision 

dynamic  frequency 

4 

E$(1l) 

single  precision 

dynamic  phase  angle 

4 

E$(12) 

single  precision 

name  of  time  history  file 

4 

E$(13) 

single  precision 

Tangential  static  load  left  side 

4 

E$(I4) 

single  precision 

static  load  right  side 

4 

E$(15) 

single  precision 

dynamic  amplitude 

4 

E$(16) 

single  precision 

dynamic  frequency 

4 

E$(17) 

single  precision 

dynamic  phase  angle 

4 

E$08) 

single  precision 

name  of  time  history  file 

4 

E$(19) 

single  precision 

Meaning  of  flag  variable: 


mu 


Distributed  transverse  static  load  (l=yes,  0=no) 
Distributed  transverse  dynamic  load  (1=yes,  0=no) 
Is  dynamic  load  harmonic?  (1=yes,  0=no) 
Distributed  tangential  static  load  (l=yes,  0=no) 
Distributed  tangential  dynamic  load  (1=yes,  0=no) 
Is  dynamic  load  harmonic?  (1=yes,  0=no) 


6  character  flag  (string  variable) 


The  above  file  structure  allows  the  user  to  apply  static  and  dynamic  loads  at  the 
same  time,  in  addition  static  distributed  loads  can  vary  linearly  (they  can  have 
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different  values  at  each  side  of  the  element).  Though  use  of  the  above  flag, 
i  element  loading  may  also  be  harmonic  or  non-harmonic. 


Displacement  History  File: 

Name:  structure  Name>.displ 
Type:  permanent,  random  access 


Info: 

Description 

Field  Field 
Length  Name 

Remarks 

displacement  of  node 

8 

n/a 

( actual  or  relative, 

) 

velocity  of  node 

8 

n/a 

( see  program  notes 

) 

acceleration  of  node 

8 

n/a 

(for  explanation. 

) 

The  purpose  of  this  file  is  to  store  the  initial  conditions  of  the  structure,  as 
defined  by  the  user  and  to  store  the  displacement  -vs-  time  history  of  the 
structure  after  solution.  Each  record  contains  the  displacement,  velocity,  and 
acceleration  of  the  appropriate  degree  of  freedom  of  the  node.  The  first  3(GN) 
(3  degrees  of  freedom  times  the  number  of  global  nodes)  records  the  initial 
conditions  of  the  structure,  t=0  (supplied  by  the  user).  The  next  3(GN)  records 
report  the  conditions  of  the  structure  at  time  step  1,  and  so  on. 

Reduction  File: 

Name:  structure  Name>.reduce 
Type:  permanent,  random  access 

Info:  Field  Field 

DEScriPtion _  Length  Name  Remarks _ 


1  Si  reduction  of  Node/DOF 

8 

n/a 

*  between  1  and  3(GN) 

reduction  of  Node/DOF 

8 

n/a 

*  between  1  and  3(GN) 

jth  reduction  of  Node/DOF 

8 

n/a 

*  between  1  and  3(GN) 

8 

n/a 

*  between  1  and  3(GN) 

The  purpose  of  the  above  file  is  to  store  a  list  of  equations  to  retain  in  the 
matrix  formulation.  The  DynFEP.create  data  file  constructs  the  above  file  with 
information  provided  by  the  user. 

To  support  the  use  of  non-harmonic  forcing  functions  or  specified  displacements 
of  nodes,  the  programs  are  capable  of  reading  a  time  history  file.  Each  record  of 
the  file  contains  a  time  and  a  magnitude  of  the  force  or  displacement.  The 
programs  will  interpolate  between  two  time  steps  if  the  required  time  is  not  on 
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file.  A  rapid  search  is  employed  by  the  program  to  find  appropriate  time,  it 
assumes  that  the  file  is  sequential  in  time.  The  first  record  in  the  file  must 
contain  the  total  number  of  time  steps  recorded  in  the  history  file. 

It  should  also  be  noted  that  the  DynFEP  program  assumes  that  information 
presented  in  this  file  is  absolute  displacement,  while  the  DynFEP.uncouple/solve 
program  assumes  that  the  information  is  accelerations. 


User  defined  Force  or  Displacement  History  File: 
Name:  <User  specified  file  name> 

Type:  permanent,  random  access 

Info:  Field  Field 

Description _  Length  Name 

Number  of  time  steps  in  file  8  n/a 

Not  used  8  n/a 


Time  8 

Displacement  or  Force  magnitude  8 


Remarks _ 

First  record  only 
First  record  only 

normal  record 


During  the  solution  of  the  problem  other  permanent  and  temporary  files  will  be 
created.  The  purpose  of  these  data  files  is  first  to  provide  storage  of  matrices 
necessary  in  the  solution  and  there  by  reduce  the  amount  of  memory  space 
required.  Secondly,  these  data  files  eliminate  the  need  to  recalculate  matrices 
to  analyze  different  loading  conditions  or  use  alternative  solutions.  A  summary 
of  these  data  files  is  presented  below. 

Permanent  data  files: 


structure 

structure 

structure 

structure 

structure 

structure 

structure 

structure 

structure 


Name>.K&F.c 

Name>.M.c 

Name>.K&F 

Name>.n 

Name>.K* 

Name>.M* 

Name>.reduce 

Name>.S 

Name>.eigen 


stiffness  and  static  forces  before  BC  3(GN)*3(GN>1 
mass  matrix  before  essential  BC  3(GN)><3(GN) 
stiffness  and  static  forces  after  BC  nx(ml) 
mass  matrix  after  essential  BC  n*n 

reduced  stiffness  and  static  forces  mx(m+1) 
reduced  mass  matrix  m^m 

listing  of  nodes  to  be  retained  mxl 

structure  mode  shapes  m*m 

eigenvalues  of  structure  n*m 
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Temporary  data  files: 

File  name _ 

structure  Name>.  Kpp 
structure  Name>.Kps 
<Structure  Name>.Ksp 
<5tructure  Name>  Kss 
<5tructure  Name>.P1 
<Structure  Name>.P2 
structure  Name>.D 


Description _ 

partitioned  stiffness  matrix 
partitioned  stiffness  matrix 
partitioned  stiffness  matrix 
partitioned  stiffness  matrix 
needed  to  calc  mode  shape  [T] 
needed  to  calc  mode  shape  [T] 
dynamic  matrix  [K]_1  [M] 


Size _ 

m*m 

m*(n-m) 

(n-m)*m 

(n-mMn-m) 

(n-m)xm 

(n-m)xm 

mxm 


Where  n  is  equal  to  the  number  of  unknown  displacements  and  m  is  equal  to  the 
number  of  modes  retained  in  the  answer.  Note  that  if  the  structure  is  not  reduced 
then  m  is  equal  to  n. 
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Conclusions 

A  set  of  micro-computer  programs  capable  of  analyzing  the  dynamic  behavior 
plane  frame  structures  has  been  developed  from  the  set  of  assumed  governing 
differential  equations.  The  system  of  programs  allow  different  approaches  to  be 
used  in  analyzing  a  dynamic  structure.  The  capability  to  use  different  approaches 
provide  a  means  of  building  confidence  by  comparing  the  results  of  the  different 
methods,  and  the  flexibilities  provided  by  the  different  approaches. 

Each  solution  approach  presented  has  unique  abilities  and  limitations.  While  the 
straight  numerical  integration  performed  by  DynFEP  has  the  ability  to  include  the 
independent  displacement  of  nodes,  it  must  process  the  formulation  with  no 
reduction  or  further  simplificatioa  The  added  capability  has  come  at  the  cost  of 
not  being  able  to  analyze  larger  structures  and  in  a  longer  computation  time.  The 
modal  analysis  approach  presented  offers  a  faster  computation  time  but 
sacrifices  the  ability  to  effectively  handle  independent  movement  of  nodes.  The 
dynamic  reduction  approach  offers  the  ability  to  do  larger  structures  with  little 
or  no  additional  computation  time,  but  at  some  sacrifice  for  accuracy. 

The  choice  of  which  solution  approach  to  use  will  normally  be  based  on  the  type 
of  problem  to  be  solved.  As  an  example,  in  Civil  Engineering  programs  such  as 
these  would  be  used  primarily  for  the  analysis  of  structure  response  to 
earthquakes.  The  dynamic  reduction  approach  presented  above  would  be  most 
useful  in  this  situation  as  it  offers  the  best  computation  time  to  size  advantage 
and  can  support  the  boundary  conditions  imposed  by  an  earthquake. 

The  programs  as  presented  here  are  capable  of  handling  about  50  nodes  when  run 
on  a  system  with  370K  of  memory  available  (assuming  8  bytes  of  memory  is 
required  for  each  double  precision  variable  used).  There  are  two  primary 
approaches  which  could  be  employed  to  increase  the  capacity  of  the  programs. 
First  the  assemblage  of  the  global  matrices  could  be  done  in  upper -banded  form. 
Since  mass  and  stiffness  matrices  are  normally  very  sparse,  this  would 
substantially  increase  the  capacity.  Secondly  the  global  mass  and  stiffness 
matrices  could  be  done  on  disk  rather  then  in  core  (also  the  application  of 
boundary  conditions,  and  the  reduction  of  the  mass  and  stiffness  matrices). 
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While  this  would  slow  down  execution  because  of  increased  I/O  activity,  the 
capability  to  handle  larger  structures  is  limited  only  by  the  amount  of  disk  space 
and  the  number  of  retained  modes.  Since  a  typical  micro-computer  system  in  a 
professional  installation  will  include  between  10  and  20  megabytes  of  hard  disk 
storage,  the  analysis  of  structures  with  several  hundred  nodes  is  seen  to  be 
reasonable. 

In  addition  to  increasing  the  capacity  of  the  programs  presented,  there  are  other 
enhancements  which  would  make  them  more  valuable  as  an  engineering  tool.  The 
programs  now  simply  grind  though  a  specified  number  of  time  steps.  It  would  be 
desirable  for  them  to  have  the  ability  to  check  selected  parameters  during  the 
processing  and  determine  for  themselves  with  to  stop  the  processing.  Such 
parameters  might  include  critical  stresses  in  selected  members,  a  maximum 
stress  in  any  member,  completion  of  a  full  cycle  of  all  forcing  functions, 
achievement  of  a  maximum  deflection,  etc. 

The  ability  to  handle  three  dimensional  plane  frames  is  a  very  natural  expansion 
to  the  program  capabilities.  Other  enhancements  could  include  an  interactive 
input  of  information  in  a  CAD/graphics  orientated  format,  and  graphic  replay  of 
structure  response. 

All  of  these  enhancements  are  within  the  current  computing  capability  of  today’s 
micro-computers.  While  the  analysis  on  truly  complex  structures  will  remain  in 
the  domain  of  mainframe  computers,  the  dynamic  analysis  of  small  structures  is 
within  the  realm  of  processing  by  micro-computer  systems.  Such  smaller 
structures  are  those  that  can  be  described  in  several  hundred  nodes  or  less,  or 
simplifications  of  more  complex  structures  which  are  being  used  for  preliminary 
design  or  investigation  prior  to  a  more  detailed  analysis  on  a  mainframe.  The 
methods  and  programs  presented  in  this  paper  form  a  corner  stone  for  building 
the  enhanced  systems  which  are  required  to  fill  this  expanding  field  of 
engineering  analysis. 
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DunFEP.menu 


Flow  Diagram  [§; 


I 


Notes: 

Program  Variables: 

GN  «  number  of  global  nodes. 

NE  -  number  of  elements, 
n  ■  number  of  unknown  displacements 
In  structure. 

m  "  number  of  retained  modes  (if  no 
reductions  then  m-n). 

Path!  >  string  defining  the  choosen 
method  of  solution. 


Valid  Menu  Selections: 

If  no  structure  data  file  is  open,  then 
valid  selections  include.  Open  File, 
and  Create  New  File, 
if  a  structure  data  file  is  opened,  then 
all  selections  are  valid. 
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'  4 - * 

'  I  DynFEP.nenu  1 

'  * - 4 


CQttON  GN,N£,DOF,n,n,PNt,Patht 
Start:  'set  up  nenu 

WINDOW  l,,<110,30)-(365,54),2:  GOSUB  BigText 

WINDOW  1:CALL  MOWETD(20 ,16) :  PRINT  ‘Dynanic  Finite  Element  Progran1; 

WINDOW  2, ,<130,70)-<350 ,306) ,4 
BUTTON  1,1, ‘Create  New  Data  Fi1e‘,(20,20)-(200,40),1 
BUTTON  2,1, ‘Open  Existing  Data  F i 1 e ' , < 20 , 45) -< 200 , 65> , 1 
BUTTON  3,0, ‘Direct  Integration  OnT y‘ ,<20 ,80)-(200 , 100) , 1 
BUTTON  4,0, ‘Nodal  Analysis*, <20, 105M200, 125),  1 
BUTTON  5,0, ‘Dynanic  Reduction", (20, 130M200, 150),  1 
BUTTON  4,1, ‘Quit*, <120, 165M2Q0, 185), 1 
IF  LB«PNt)>Q  THEN  GOSUB  Look.at.File 
60SUB  NornalText 

loop: 

CALL  M0UET0< 10 ,200) :  PRINT  ‘Problen  nane  =  ‘;PN* 

CALL  NOVETO<10,212):  PRINT  USING  ‘Global  Nodes  =««’ ;GN 
CALL  NOWETOd 0,224):  PRINT  USING  ‘Elenents  =«‘;NE 
WHILE  DIALOG<0)<)1:  UEND  'wait  for  the  user  to  do  sonething 

ON  OlALOG(l)  GOTO  Create, Existing, FEN, Modal , Reduce, Quit  'branch  according  to  nenu  selection 

Create:  Patht=*Create  New  Data  File*:  GOSUB  Change .Window 
Ft=*Basic  Disk  l:DynFEP. create  data  file* 

CHAIN  Ft:  BID 

Existing: 

PNt=FILESt<l ,‘DYNA‘)  '  get  Problen  Nane  f ran  user 
IF  PNt<>‘‘  THEN  GOSUB  Look.at.File 
GOTO  loop 

FEN:  Patht=‘Direct  Integration  Only1:  GOSUB  Change .Window 

IF  n<0  THEN  Ft=’Basic  Disk  l:DynFEP.nass/stiffness‘  ELSE  Ft=‘Basic  Disk  l:DynFEP‘  'no  need  to  reassenble  gobal  natric 
es. 

CHAIN  Ft:  END 

Nodal:  Patht=‘Modal  Analysis*:  GOSUB  Change .Window 

IF  n<0  THEN  Ft=*Basic  Disk  l:DynFEP.nass/stiffness‘  ELSE  Ft=*Basic  Disk  1  :DynFEP. essential  BC*  'no  need  to  reassenble 
gobal  natrices. 

CHAIN  Ft:  END 

Reduce:  Patht=‘Dynanic  Reduction’:  GOSUB  Change .Window 

IF  n<0  THEN  Ft=‘Basic  Disk  1 iDynFEP.nass/stiffness*  ELSE  Ft=‘Basic  Disk  1 :DynFEP. essential  BC‘  'no  need  to  reassenble 
gobal  natrices. 

ClttlN  Ft:  END 

Quit:  Patht=‘0utput  Window*:  GOSUB  Change .Window 
BID 


Subroutines  Below 
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Look.At.File:  D0F=3  '  Find  status  oF  processing 

0PB4  PNi  FOR  INPUT  AS  (1:  1NPUTI1 tGN,NE,nfn:  CLOSE  II 
FOR  i=3  TO  5:  BUTTON  i,l:  NEXT  i  'activate  buttons 
RETURN 


Change .Window:  UINDOU  CLOSE  1:  WINDOW  CLOSE  2:  CLOSE 
WINDOW  I ,Path<,<5,40)-(265,298),l :  RETURN 

BigText :CALL  TEXTF0NT<0) :CALL  TEXTSIZE(12):RETURN  '  Chicago 
LittleText:CALL  TEXTF0NT(1):CALL  TEXTSI2E(9) :RETURN  '  Geneva 
NomalText  :CALL  TEXTFONT <  1 ) : CALL  TEXTSI 2E<  10): RETURN  '  Geneva 
FornatedTextiCALL  TEXTF0NT(4):CALL  TEXTSIZE(9):RETURN  '  Monaco 


i 

I 


i 


1 
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Craata  nod*,  alvnant.  and  Info 
data  run 


j  WMdrad  Man  flag 


I  Mad  raductlon  data  I 

.  ..r 

I  Star*  reduction  data 


Notes: 

Program  Variables: 

6N  -  number  of  global  nodes  .s 
NE  *  number  of  elements, 
n  -  number  of  undnown  displacementss 
in  structure. 

m  -  number  of  retained  modes  (if  no 
reduction  then  m=n) 

Path$  -  string  variable  indicating  the 
choosen  method  of  solution. 

The  program  assumes  that  a  data  Text 
file  has  been  prepared.  Format  of  the 
file  is  that  of  a  BASIC  DATA  statement 
All  data  shown  in  the  discription  of 
the  data  Hie  structure  must  be  included. 
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'  I  DynFEP. create  data  tilt  I 


C0N10N  GN,NE, DOF, n,ii,PN4, Path* 

WINDOU  2 , , < 1 1 0 , 250 ) -( 380 , 300 ) , 2 :  GOSUB  BigText:  U1NDQU  2 
PRINT  'It  you  have  not  created  a  text  tile  ot* 

PRINT  'Basic  DATA  statenents  tor  this  progran,' 

PRINT  'press  the  'Cancel'  button ! "CHW<7) ; 

PN$=F1LESK1,'TEXT')  'get  data  tile  nane  trai  user 
WINDOU  CLOSE  2:  IF  PN$="  THEN  STOP  'user  hit  cancel  button 
PRINT  'Type  'GOTO  Start'  to  continue. *CHR«7) 

MERGE  PN$  'load  user  created  data  tile  &  execute  with  it 

Start:  DIM  N$(17),E«19):ZX=-1 
GOSUB  FomatedText 
READ  PN$ 

F*=PNK'. nodes':  OPEN  Ft  AS  II  LEN=92 

FIELD! 1 ,12  AS  FIgl«,  4  AS  NK1),  4  AS  NK2),  4  AS  NK3),  4  AS  NK4),  4  AS  Nf<5),  4  AS  NKd),  8  AS  NK7),  4  AS  NK8),  4  A 

S  N$(9) ,  4  AS  NK10),  4  AS  N«ll),  8  AS  N$(12),  4  AS  NK13),  4  AS  NK14),  4  AS  Ni(lS),  4  AS  NK16),  8  AS  NK17) 

F<=PNK'.e1enents*:  OPEN  F$  AS  12  LEN=94 

FIELDN2 ,6  AS  F1g»,2  AS  Lt$,2  AS  Rt*,4  AS  ES(1) ,4  AS  EK2),4  AS  EK3),4  AS  EK4),4  AS  EK5>,4  AS  EK4),4  AS  EK?),4  AS  E 

1(8), 4  AS  EK9),4  AS  E«10),4  AS  EK11),4  AS  EK12),8  AS  EK13),4  AS  EK14),4  AS  EK15),4  AS  E$(li),4  AS  E$(17),4  AS  EK 

18) ,8  AS  E«19) 

F«=P9»:  OPEN  F$  FOR  OUTPUT  AS  15 

READ  NunNodesX 
FOR  i=l  TO  NunNodesX 
READ  a*:  LSET  F1gl4=at 
FOR  j=l  TO  17 

IF  js7  OR  j=12  OR  j=17  THEN  READ  a«:  LSET  NKj)=al  ELSE  READ  a:  LSET  NKj)=HKSKa) 

NEXT  j 
PUTI1 ,  i 
NEXT  i 

READ  NunElementsX 
FOR  i=l  TO  NunElenentsX 

READ  al,LtX,RtX:  LSET  Flg2$=al:  LSET  Lt*=MKlKLtX):  LSET  R«=MKH(RtX) 

FOR  j=l  TO  19 

IF  j=13  OR  j=19  THEN  READ  af:  LSET  EKj)=a$  ELSE  READ  a:  LSET  EKj)=MKSKa) 

NEXT  j 
PUTI2, i 
NEXT  i 

'  HHHHHHHIHHHHHHHIHHH  debug 

'TRON 

'HIHHHIHIlllliHIHIilHIIIlHH 

DIM  UI(NunNodesX*3,3)  '  initial  conditions 
FOR  i-1  TO  NunNodesX*3 

READ  UI<i,l),UI(it2),UI(i,3) 

NEXT  i:  n>3:  n=NunNodesX*3 

CALL  Display .Matrix(n,n3, UK), 'Initial  Conditions') 

CALL  store  .Matrix(n,n3, UK)  ,PNK'.  ini  tial',n3) 

'  IIHIIHHHHHIliHHHHHHIHH  debug 
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'TROFF 

'WHWHHHHHHiHHHmHHH 

INPUT  ’Reduction  into  in  this  data  <y/n)*jal 

bfO:  IF  aRO'y*  THEN  6OT0  finish. up 

READ  n:  QPB4  PNR+*. reduce*  AS  14  LB4=8:  FIELD«4,8  AS  aat 

FOR  i*l  TO  m:  READ  r:  LSET  aa<=HKM(r):  PUT«4,i:  NEXT  i 

finish. up; 

INPUT  ‘Have  global  matrices  been  assembled  <y/n)*;at 
IF  a«OV  THEN  n=-l  ELSE  n=0 
WRITERS  ,NunNodesX  ,NumEl  emen  t  sX ,  n  ,m 
GN=NunNodesX:  NE=NunElenentsX:  00F=3 
CLOSE:  NAME  PNS  AS  PN4,'DYNA' 
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DunFEP-mma/tt  1  rfnns 

_ Flow  Diagrem _ 


Notes: 

Program  Variables: 

6N  -  number  of  global  nodes. 

NE  -  number  of  elements. 

n  -  number  of  unknowns  in  the  structure. 

m  -  number  of  retained  modes  (if  no  reduction 
then  m-n). 

Path!  •  string  variable  indicating  the  choosen 
method  of  solution. 

[K  J  provides  storage  for  the  global  stiffness 
matrix,  and  the  static  force  matrix. 

The  stiffness  matrix  is  a  square  matrix 
with  dimensions  equal  to  GN  times  3. 

The  static  force  matrix  is  stored  with 
the  stiffness  in  an  addition  column. 

[Ml  provides  storage  for  the  global  mass 
matrix.  It  is  a  square  matrix  with 
dimensions  equal  to  GN  times  3. 
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'  I  DynFEP.mass/stiffness  I 

'  4 - 4 

CtttIGN  GN, NE, DOF, n,m,PN$, Path*:  GOTO  Start 


Subroutines  below 


AssembleForceNat: 

FOR  i=l  TO  GN:  6ETI1 , i  '  read  specified  nodal  loads 
FOR  j=0  TO  DOF-1 
i ndex= < i - 1 ) »OOF+ j ♦ 1 : k= 3+ 5» j 
' — NodeStatForces: 

Sload=0:  A=l*j*4 

IF  THEN  Sload=CVS<N*(k))  '  static  load 

K#< index, n+l)=K#( index ,n+l)+S1oad  '  add  in  force;  positive  to  right,  upward,  clockwise 

NEXT  j,i:  RETURN 

El ementNatr  ixAssembl er : 

FOR  ELENENT=1  TO  NE 
60SUB  Bui IdElenentliatr ices 
1R=(N1X-1)»D0F:IC=(N2X-1)*D0F 
' — Assen.K.M.Stat. Elen. Forces: 

FOR  i=l  TO  OOF 

K*(IR+i,n+l)=KI(IR*i,n+l)4A»(i,7)  '  elenent  forces  stored  in  column  »n+l 
K#< IC* i ,n+l)=Kf< IC+ i ,n+ 1 HA#< i +D0F ,7) 

FOR  j=l  TO  OOF 

KI<IR+i,lR+j)=K*<IR+i,IR*j)+AR(i,j)  '  assemble  stiffness  matrix 
K*(IR+i ,IC+j)=K»(IR+i ,IC+j)+A»(i ,j*D0F) 

KICICe  j  f  IR-t  j)=Ki(ICe  i  ,IRe  j  )4Ai<  i +D0F ,  j ) 

KI(ICei,lC*j)=KI(IC+i  ,IC+j)+AII(itDOF,j+DOF) 

NI(IR+i,IR+j)=MR(IR+i,IR*j)*BI(i,j)  '  assemble  mass  matrix 
M«<IR+i  ,IC*j)=MKIR*i  fIC*j)*B»<i  ,j+D0F) 

Mi<IC*  i ,  IR*  j)=til<IC*  i  ,IRe  j  )*BS<  i*D0F,  j> 

HKIC+i  ,IC+j)st1l(IC+i  ,IC*j)+BI(i+D0F,j+D0F) 

NEXT  j,i 

NEXT  element 

Bui  ldEleiaentHatr  ices: 

GETI2, ELENENT :NlZ=CVI<Ltt):N2X=CWI(Rti)  '  get  left  and  right  global  node  ITs 
GETI1 ,N1X :X1=CVS(N$( 1 ) ) :Y1=CVS(N*<  2) )  '  get  left  side  coord's 
GETI 1 ,N2X :X2=CVS(N$( 1 ) ) :Y2=CVS<N$ ( 2) )  '  get  right  side  coord's 
L=SQR«Y1-Y2)*2+<X1-X2)*2)  '  find  element  length 
GOSUB  Elem.K.M.Stat .Forces 

IF  Xl-X2=0  THEN  angle=SGN<Yl-Y2)*Pil/2  ELSE  angle=2*Pi»-ATN((Yl-Y2)/(Xl-X2)) 

IF  angle>2*Pi«e .003  OR  angle<2«Pi»-.003  THEN  60SUB  Transform 
RETURN 

Elem.K.N.Stat .Forces: 

I1=CWS(E*(1)):I2=CVS(E«4))  '  moments  of  inertia 
A1=CWS< E*( 2) )  sA2=CVS< E*<  5) )  '  areas 
ml=CVS(E*<3)):m2=CVS<E*<4)>  '  mass/length 
E=CVS<E<<7))  '  elastic  modulus 
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FOR  i=l  TO  4:F0R  j=l  TO  7:AI(i ,j)=0:NEXT  j,i  '  initialize/build  element  stiffnesses 
Al< 1 , 1 )*E*(A1 +A2)/(2*L) :AI<4 ,4)=A»( 1 , 1 ) :A«( 1 ,  4)=-AI< 1,1) 

AI<2,2)=E*4*dl+I2)A43:AI<5,5)=AK2,2):AK2,5)=-AI<2,2) 

AI<2,3)=-E*<4«I1+2*I2)A42:AK3,5)=-AI<2,3) 

AK  2 , 4)=-E*<2*l  1 +4*1 2)  A4  2  :AI<  5 , 4)=-AI<  2 ,4) 

AK3,3)SE*<3*I1+I2)A 

AI<3,4)=E*dl+I2)A 

AI(4,4)=E*dl+3*I2)A 

FOR  i=2  TO  4:F0R  j=l  TO  i-1  :AI< i ,j)=Ai< j , i ) sNEXT  j,i  '  symetrize  st i-f-fness 

FOR  i=l  TO  4:F0R  j=l  TO  4:BI<i,j)=0:NEXT  j,i  '  initialize/build  element  matrix 

Bid  ,l)*70*L*<3*ml+m2)/840 

Bid  ,4)=70*L*<ml+m2)/840 

B» ( 2 , 2)=24*L*< 1 0<ml +3*m2)/840 

B»( 2 ,3)=-2*L‘2*< 15*ml +7*m2)/840 

BK2,5)'54*L«ml+m2)/840 

BK 2 , 4)=2*L4  2*<  7*ml +4*m2)/840 

BK  3 , 3)=L 4  3*  ( 5*ml + 3*m2 )/840 

Bl< 3 , 5)=-2*L4  2*<i*ml +7*m2)/840 

8K3,4)=-3*L43*<ml+m2>/840 

B»( 4 , 4)=70*L*(ml +3*m2)/840 

B»<5,5)=24*L*(3*ml+10*m2)/840 

BK5 , 4)=2*L4  2*<  7*ml  ♦  1 5*m2)/840 

BI(6,4)=L43*(3*ml+5*m2)/840 

FOR  i=2  TO  4:F0R  jsl  TO  i-l:BKi,j)=BKj,i):NEXT  j,i  '  symetrize  mass  matrix 
DSloadl=0:DSload2=0:TSloadl=0:TSload2=0  '  Find  element  loading 

IF  MIDKFlg2*,l,l)=M*  THEN  DSloadl=CVS<EK8)) :0S1oad2=CVS<EK9»  '  distributed  static  load 
IF  IUW(Flg2$,4,l)=,l*  THEN  TSloadl=CUS(Ef<14):TSload2aCVS(EI<15)  '  tangential  static  load 
Aid ,7)=L*(20*TS1oadl+10*TSload2)/40  '  positive  to  the  right 

AK2,7)=L*<-15*DSloadl+45*DSload2)/40  '  positive  upward 

AK3,7)-*L42*(3*DSloadl+2*DSload2)/40  '  positive  clockwise 

AK4,7)=L*d0*TSloadl+20*TSload2)/40  '  positive  to  the  right 

AI(5,7)=L*<9*DSloadl+21*DSload2)/40  '  positive  upward 

AK4,7)=L42*(2*DSloadl+3*DSload2)/40  '  positive  clockwise 

RETURN 

Transform:  'Subroutine  to  transform  Stiffness,  Mass,  and  Force  element  matrices 
GOSUB  Bu i 1 dTransf ormat i onMat 

FOR  i=l  TO  4:F0R  j=l  TO  4:CKi,j)=0:FOR  k=l  TO  4:CKi,j)=CKi,j)+TI(k,i)*AI<k,j):NEXT  k,j,i  'transposelTWKel 

FOR  i=l  TO  4:F0R  j=l  TO  4:AKi  ,j)=Q:F0R  k=l  TO  4:AKi,j)=AKi,j)+CKi,k)*TI(k,j):NEXT  k,j,i  '[transposelTWKeimT] 

FOR  i=l  TO  4:CKi,l)=0:F0R  k=l  TO  4 :C*< i , 1 )=CI< i ,  1 )+TI<k , i )*AI(k ,7) :NEXT  k,i:F0R  i=l  TO  4:AKi,7)=CKi,l):NEXT  'trantTl 
»<Fe) 

FOR  i=l  TO  4: FOR  j=l  TO  4:CKi,j)=0:F0R  k=l  TO  4:CKi,j)=CKi,j)+TKk,i)*BKk,j):NEXT  k,j,i  'transposelTHIMel 

FOR  i=l  TO  4:F0R  j=l  TO  4:BKi  ,j)=0:FOR  k=l  TO  4:BKi,j)=BKi,j)+CKi,k)*TKk,j):NEXT  k,j,i  'ttransposem*[Me]]*[T] 

RETURN 


BuildTransformationMat:  '  build  IT) 

FOR  i=l  TO  4:F0R  j=l  TO  4:T»(i ,j)=0 :NEXT  j,i  '  initialize 
IF  angle  MOO  Pil/2  THEN  Tld,l)=COS(angle)  ELSE  TI<1 ,1)=0 
Tl(  4,4)=TI<  1 , 1 )  :TI<  2 , 2)=TI(1 , 1 )  :TI<  5 , 5)=Tld ,  1 ) 

IF  angle  MOD  Pil  THEN  Tld ,2)=-SJN<angle)  ELSE  7ld,2)=0 
TI(4,5)=TI<l,2):TI(2,l)=-Tld,2):TI(5,4)=TK2,l) 
TK3,3)=1:7I(4,4)=1 
RETURN 
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Start:  n=GN*D0F:  al=l:  a3=3:  CRt=CHRt<13):  PilM*ATN(l) 

DIH  KI(n,n+l),NI(n,n) ,AI<6,7) .81(6,6) ,0(6,6) ,T«<6,6) ,Nt< 17) ,Et(19) 

Ft=PNt4*. nodes*: OPEN  Ft  AS  »1  LEN=92 

FIELDtl ,12  AS  F1glt,4  AS  Nt(l),4  AS  Nt(2),4  AS  Nt(3),4  AS  Nt<4),4  AS  Nt(5),4  AS  N$<6),8  AS  Nt(7),4  AS  Nf<8),4  AS  Nt(9),4 
AS  Nt(10),4  AS  Nt(ll),8  AS  Nt<12),4  AS  Nt<13),4  AS  Nt<14),4  AS  Nt(15),4  AS  Nt<16),8  AS  N$<17) 

Ft=PNt+ '.elements': OPEN  Ft  AS  12  LEN*94 

F1ELD#2,6  AS  Flg2t,2  AS  Ltt,2  AS  Rtt,4  AS  Et(l),4  AS  Et<2),4  AS  Et<3),4  AS  Et<4),4  AS  Et<5),4  AS  Et(6),4  AS  Et(7),4  AS  E 
1(8) ,4  AS  Et(9>,4  AS  Et(10),4  AS  Et(ll),4  AS  Et(12),8  AS  Et(13),4  AS  Et(14),4  AS  Et(15),4  AS  E«(16),4  AS  Et(17),4  AS  Et( 
18), 8  AS  Et<19) 

Build. Global  .Na trices: 

60SUB  ElementNatrixAssembler 
ERASE  At,Bt,CI,T» 

Build. Static. Force.Nat: 

DIN  U0i(n,3) ,Ul»(n,3) ,Qt(n) 

GOSUB  AssenbleForceNat 

CALL  Store.Natrix<n,n+l,K»(),PNt+'.K4F.c*,a3)  '  store  stiffness  and  force  matrices 
CALL  Store.Natrix<n,n,H»0,PNt+\N.c*Ia3)  ’  store  stiffness  and  force  matrices 

'iHIHHHlHHHIIIHMHHIiHIHHHI  '  debug 

CALL  Display .Matrix(n,n+1,K#(), 'Stiffness') 

CALL  Display .Matrixtn.n.MIO, ’Mass') 

'HiHHIHHlIlHHHHHHIHHiHIHH 

n=fl:  OPEN  PNt  FOR  OUTPUT  AS  13:  WRITE*3,GN,NE,n,m:  CL0SEI3:  NONE  PNt  AS  PNt,*DYNA* 

CLOSE 

IF  Patht=*Direct  Integration  Only*  THEN  CMIN  ‘Basic  Disk  1  :DynFEP*  ELSE  CHAIN  ’Basic  Disk  1 :DynFEP. essential  BC* 
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Notes: 

Program  Variables: 

GN  -  number  of  global  nodes. 

NE  -  number  of  elements. 

n  -  number  of  unknowns  in  the  structure. 

m  -  number  of  retained  modes  (if  no  reduction 
then  m=n). 

Path!  -  string  variable  Indicating  the  choosen 
method  of  solution. 

[K]  provides  storage  for  the  global  stiffness 
matrix,  and  the  static  force  matrix. 

The  stiffness  matrix  is  a  square  matrix 
with  dimensions  equal  to  GN  times  3. 

The  static  force  matrix  is  stored  with 
the  stirrness  in  an  addition  column. 

(Ml  provides  storage  for  the  global  mass 
matrix.  It  is  a  square  matrix  with 
dimensions  equal  to  GN  times  3. 

(BC)  provides  storage  for  the  boundary 
condition  index,  a  column  of  1/0s. 

If  0  then  displacement  has  been  specified. 
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'  * - 4 

'  I  DynFEP.essential  BC  ! 

'  4 - 4 


C0N1QN  GN,NE, DOF, n,i»,PN$, Path* 

'  Subroutines  Below 


Suitch: 

FOR  j=l  TO  n:  SUAP  K«< i ,j) ,K»<K, j) :  SUAP  MUti  ,j),MKk,j):  NEXT  j 
FOR  j=l  TO  n:  SUAP  K#<j,i),K#<j,k):  SUAP  M#(j,i)pM#(j,k):  NEXT  j 
RETURN 


Starts  n=GN*DOF:  al-1 :  a2=2 

DIM  8Ci<n,l) ,K#<n,n4l) fMiCn,n) ,N»<17) 

F$=PN$4*. nodes* sOPEN  F$  AS  HI  LEN=92 

FlELDtll ,12  AS  Flgl$,  4  AS  N$(l),  4  AS  N*<2),  4  AS  N*<3),  4  AS  N$<4),  4  AS  N*<5),  4  AS  N$<6),  8  AS  N*<7),  4  AS  N$<8),  4  A 
S  N$(9),  4  AS  N$(10),  4  AS  N$(ll),  8  AS  N$(12),  4  AS  N$(13),  4  AS  N*(14),  4  AS  N$<15),  4  AS  Ni(lA),  8  AS  NS<17) 

CALL  Retrieve .Matrix(n,n4l,K«()lPN$4M(&F.c*la2) 

CALL  Retrieve.Matrix(n,n,MIK),PN$4*.M.c*,a2) 

Detine. Essential .BC: 

FOR  i=l  TO  GN:  GETll.i 
FOR  j=l  TO  DOF 

index=(i-l)*DOF4j:  k=l44»<j-l) 

BCIK  index,  1)=UAL(M1D$(F1  gif, k,D) 

NEXT  j,i 

'  HHHHHHIHHHHiiHHlHiHHHH  debug 

CALL  Display .Matrix<n,al,BdK),*8.  C.*) 

CALL  Di splay .Matrix(n,n4i ,KM<) , * Stiffness*) 

CALL  D i sp 1  ay .Matr i x ( n , n ,M#< ) , 'Mass* ) 

CALL  Store. Matr ix<n,al, BCIK  ),PN*4*.BC*,a2) 

Apply. Essential ,BC:  k=0 
FOR  i=l  TO  n 

IF  BCKi,l)=l  THEN  k=k4l:  IF  i Ok  THEN  GOSUB  Switch 
NEXT  i 

k=0:  FOR  i=l  TO  n:  k=k4BCi< i ,1):  NEXT  i:  n=k 
FOR  i=l  TO  n:  SUAP  K»<i,n4l),K!Ki)GN*D0F4l):  NEXT  i 
'  IHIHHiliHHlillHHIiiHHHiHHH  debug 

CALL  Display .Matrix<n,n4l,K»<), ’Stiffness*) 

CALL  D  i  sp  lay  .Matr  i  x<  n  ,n  fHtl< ) ,  ‘Mass* ) 

'HHHHIHHHiHHHHHHHHHimi 

CALL  Store .Matr ix<n,n4l ,KI(),PN$4*,K4F* ,a2) 

CALL  Store. Matr ix(n,n,MIK),PN$4*.M*,a2) 

IF  Path*=*Modal  Analysis*  THEN  npn  'if  not  n  has  been  set  by  DynFEP. create 
OPEN  PN*  FOR  OUTPUT  AS  12:  URlTE#2,GN,NE,n,m:  CL0SE42:  WME  PNi  AS  PN$,*DYNA*:  CLOSE 

IF  Path*=*Modal  Analysis*  THEN  OWN  ‘Basic  Disk  1 : DynFEP. eigen  solver*  ELSE  CMIN  ’Basic  Disk  1  sDynFEP.reduce* 

END 


-SUB  Retrieve.Matrix(r,c,AIO,F*,k)  •' 


.*« 
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DunFEP.reduce 

Flow  Diagram 


Notes: 


Program  Variables: 

GN  =  number  of  global  nodes . 

N£  -  number  of  elements. 

n  =  number  of  unknowns  in  the  structure. 

m  *  number  of  retained  modes  (if  no  reduction 
then  m=n). 

Path$  -  string  variable  indicating  the  choosen 
method  of  solution. 

IK]  provides  storage  for  the  global  stiffness 
matrix,  and  the  static  force  matrix. 

The  stiffness  matrix  is  a  square  matrix 
with  dimensions  equal  to  GN  times  3. 

The  static  force  matrix  is  stored  with 
the  stiffness  in  an  addition  column. 

[li]  provides  storage  for  the  global  mass 
matrix.  It  is  a  square  matrix  with 
dimensions  equal  to  GN  times  3. 

(R)  provided  storage  for  the  reduce  index 
It  is  a  list  of  equations  to  be  retained. 

[PI  ]  and  IP2]  are  calculated  and  stored  to 
for  use  by  eigen  solver  in  transforming 
eigen  vectors. 


Available  Sub-Programs: 

Display  .Matrix 
Store  .Matrix 
Retrieve  .Matrix 
Mat. time  .Mat 
Mat  .pi  us  .Mat 
lnvert.Matrix 
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'  ♦ - + 

'  i  OynFEP. reduce  I 

'  4 - 4 


COftlON  GN , NE, DOF, n,n,PN$, Path* 

Starts  al=l:  p=n:  s=n-n:  nll=l 

DIM  Kh(n,n+1) ,Ri(n,l)|Kpptt(p,p) ,KssKs,s) ,KspKs,p) , Kps#< p , s) :  IF  p>s  THEN  d=p  ELSE  t*s 
DIM  Tl#<d,d) ,T2S<d,d)  'tenporary  storage 
CALL  Retr ieve. Hatr ix<n,n*l ,Ki<) ,PN*+' .K4F’ ,al) 

CALL  Retr ieve.Natrix<n,al,RK),PN*4a. reduce1, al) 

'  liHiHIHHilMHHHHHHHiHHHiHH  debug 

CALL  Display .Matr ix<n,al ,RiO , 'Equations  to  be  Retained') 

'HHHHKHHHiiiHHHiHHIHHHHiH 

FOR  i=l  TO  n  '  nove  equations  to  be  retained  to  the  top 
IF  iORKi  ,1)  THEN  FOR  j=l  TO  n:  SWAP  KKi,j),KKRKi,l),j):  NEXT  j  'swap  row 
IF  iORKi, 1)  THEN  FOR  k=l  TO  n:  SWAP  K«<K , i ) ,K»< K ,Rtt< i ,1)):  NEXT  k  'swap  colunn 
NEXT  i 

FOR  i=l  TO  p:  FOR  j=l  TO  p:  KppKi,j)=KKi,j):  NEXT  j  'build  partitioned  natrices 
FOR  k=p4l  TO  n:  KpsKi  ,k-p)=KKi  ,k):  NEXT  k,i 
FOR  i=p4l  TO  ns  FOR  j=p*l  TO  n:  KssKi-p,j-p)=KKi,j>!  NEXT  j 
FOR  k=l  TO  p:  Kspt<i-p,k)=K#<i,k):  NEXT  k,i 

'HiiiiiiHiHHiHiHHHlHHtHHHHiHHiHliHiHliHiHHiiHiiHH  debug 
CALL  Display. Matrix<p, p, KppIO, aKppa):  CALL  Display .Natrix(p, s, KpsK) , 'Kps') 

CALL  OispIay.Natrix<s,p,KspK),*Kspa):  CALL  Display .Hatrix<s,s,KssK),aKssa) 

'HHHHHHmHHHHHIHiHttiHilHHHHHHHHiHHHHHiHiHi 

CALL  lnvert.Natrix<s,Kssl())  'Find  [Kssl  inverse  then  save  partitioned  natrices 
CALL  Store . Hatr  ix(p, p, KppK )  ,PN*+a .Kpp* ,al):  CALL  Store .Natrix(p,s,Kps«0,PN*4*. Kps’, al) 
CALL  Store. Hatr ix(s,p,Kspl<),PN$4*.Ksp*,al)i  CALL  Store.Natrix<s,s,KssK),PN*4\Kssa,al) 

CALL  Mat.tines.Nat<5,p,s,KssK),KspK),TlK))  'Find  the  reduced  stiFFness  natrix 
CALL  Nat.tines.Nat<p,p,s,Kpstt<),Tll<),T2«)) 

CALL  Nat. plus. Nat(p,p, nil, KppK), -nil, T2IO) 

'  HHHIHHHlHHHHiHHilHHHIHHH  debug 

CALL  Display .Natrix(p,p,Kppl(), ‘Reduced  StiFFness*) 

'iiHHIlHHHlHHHHHIHIHIHHHHH 

CALL  Store.Natrix(p,p, KppIO, PN*4*.K*',al)  'store  the  reduced  stiFFness  natrix 
ERASE  Kl , Kppl , Kspl , Kpsl , Kssl 

DIN  Nl(n,n) ,Nppl<p,p) ,Nssl(s,s) ,Nspl<s,p) ,Mpsl(p,s) 

CALL  Retr i eve .Natr i x(n ,n ,NIO  ,PNKa  ,N* , al ) 

FOR  i=l  TO  n  '  nove  equations  to  be  retained  to  the  top 
IF  iORKi, I)  THEN  FOR  j=l  TO  n:  SWAP  Ml< i , j) ,MKR«< i ,1) , j) :  NEXT  j  'swap  row 
IF  iORKi, 1)  THEN  FOR  k=l  TO  n:  SWAP  NI<k,i),MKk,RKi ,1)):  NEXT  k  'swap  colunn 
NEXT  i 

FOR  i=l  TO  p:  FOR  j=l  TO  p:  NppKi  ,j)=NKi  ,j):  NEXT  j  'build  partitioned  natrices 
FOR  k=p4l  TO  n:  NpsKi,k-p)=MKifK):  NEXT  k,i 
FOR  i*pO  TO  ns  FOR  j=p«l  TO  n:  NssKi-p,j-p)=NKi,j):  NEXT  j 
FOR  k*l  TO  p:  Mspl<  i -p ,K)=til<  i  tk ) :  NEXT  k,i 
ERASE  Nl 
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debug 

CALL  Display .liatr i x<p ,Hppf 0 taMppa) :  CALL  Display .Hatrixtp,s,MpsiO,,Mps*) 

CALL  Display .Matrix<s,p(Hsp<<) ,aMsp'> :  CALL  Display .Matrixts.SjMsslOj'Mss*) 

'HIHHIHillHimiHliHHliiHHimHHHimiHHiHHIHiHiWIH 

DIM  Kppl<p,p),Kssl<s,s),Ksplts,p),Kpsi<p,s)  'reload  the  partitioned  stiFFness  matrices 
CALL  Retrieve .Matr ix<p ,p ,Kppfl< ) ,PN*+* .Kpp‘ , al ) :  CALL  Retrieve .Matr ix(p ,sfKpsVO ,PNf+* . Kps’ ,al) 

CALL  Retrieve.Matrixts,p,Ksplt),PNi+,.Ksp,,al):  CALL  Retrieve .Matrixts^KsslOtPNi+’.Kss*  ,al)  'recall  stored  CKssl  inve 
rse 

'  Find  the  reduced  mass  matrix  and  CPI 1  and  CP23  For  use  in  Finding  CTJ  by  OynFEP. eigen  solver 
FOR  i=l  TO  d:  FOR  j=l  TO  ds  Tllti,j)=0:  T2»< i ,j)=0:  NEXT  j,i  ' init  Tl»  and  T2I 
CALL  Nat. times. Hat (s,p,s,Ksslt),Ksplt), Tilt) )  '  [Til  =  CKsslinvEKspl  =  EPU 
CALL  Store.Natrix(s,p,TllO,PN$+,.Pr,al)  'used  by  DynFEP. eigen  solver 
CALL  Mat. times.Mat[p,p,s,Mpslt),Tll<),T2IO)  '  IT2)  =  [MpsHKss] invLKsp] 

CALL  Mat. plus.Mattp.p, nil, MppIO, -nil, T2IO)  '  EMppl  =  EMppl  -  EMpslEKsslinvEKsp] 

FOR  i=l  TO  d:  FOR  j=l  TO  d:  T2lti,j>=0:  NEXT  j,i  'init  T2I 

CALL  Mat. times. Matts, p,s,Mssl(),Tll(),T2l())  '  ET21  =  EHsslEKssl invEKsp] 

CALL  Mat .plus.Mat<s,p ,nll,Mspl< ) ,-nll ,T2I< ) )  '  (Mspl  *  CMspI  ♦  EHsslEKssl invEKsp] 

FOR  i=l  TO  d:  FOR  j=l  TO  d:  Tll<i,j)=0:  NEXT  j,i  'init  Til 

CALL  Mat. times.Mat<s,p,s,Kssl<),T2IO,TllO)  '  [Til  =  EKsslinvEMsslEKssl invEKsp] 

FOR  i=l  TO  d:  FOR  j=l  TO  d:  T2Ki,j)=0:  NEXT  j,i  'init  T2I 

CALL  Mat . t imes .Nat <  s ,p , s , Ksslt ) ,Mspl< ) ,T2I< ) )  '  ET21  =  EKsslinvEHspl 

CALL  Mat. plus. Matts, p, nil, T1IO  ,-nll, T2IO)  '  [Til  =  -EKsslinvEHspl  +  EKsslinvEMsslEKssl  invEKsp]  =  EP21 
CALL  Store .Natrixts,p, Tilt), PN$4‘.P2‘,al)  'used  by  DynFEP. eigen  solver 
FOR  i=l  TO  d:  FOR  j=l  TO  d:  Tll<i,j)=0:  T2lti,j)=0:  NEXT  j,i  'init  Til  and  T2I 
CALL  Mat. times. Matts, p,s, Ksslt ),Msplt), Tilt))  '  ET11  =  [KsslinvtEMspl  ♦  EMsslEKssJinvEKspl) 

CALL  Mat. times.Mattp,p,s, Kpslt), Tilt), T2IO)  '  ET21  =  EKpslEKsslinvtEMspl  ♦  EMsslEKsslinvEKspl) 

CALL  Mat.plus.Mat(p,p, nil, MppIO, -nil, T2lt))  '  EMppl  =  EMppl  -  EMpslEKsslinvEKsp]  -  EKpslEKsslinvtEMspl  ♦  [MsslEKsslin 
vEKspl) 

'  HHHHHHiHlHHHHHHIIIHHIHHH  debug 

CALL  Display .Matrixtp,p, MppIO, 'Reduced  Mass') 

'IHliilimilllillHIHIIlHHHHHHHiH 

CALL  Store. Matr ix(p,p, MppIO, PN$+'.M»,,al)  'store  the  reduced  mass  matrix 
ERASE  Mppl,Hpsl,Mspl,Mssl 

CLOSE:  KILL  PNI+'.Kpp’:  KILL  PNf+'.Kps*:  KILL  PNIe'.Ksp':  KILL  PNH'.Kss'  'distroy  temporary  Files 
CHAIN  ‘Basic  Disk  1: DynFEP. eigen  solver* 

END 
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Flow  Diaaram 


UM  mm!  (tirriMM  matrix 
UO 


IMP  (Itrrxaaa  matrix 
HO 


Kwart  ttlffmM  matrix 
ataratnOO 


Fine  ID)  ■  Odlnaaraa  ■  Ml 


ClaarlK) 


Load  IPI1  and  IP2) 


Find  impravad  mada  ■ 

waa.PKnau 


Find  alganaalua  and 
normallM  mada  ahapa 
<naa)  a  tma)/aipan*3 


nancata  good  mada 


9tncun 


Find  tranafarmatian  matrix 
CT| » iru  ♦  aiganytrai 


lro2) .  rrxran 


St  ora  (hmi  aa  column  I 

mada  aitapa  matrix 


Siara  aifiwalua  in 
aifanaalaa  matrix 


Program  Variables: 

GN  *  number  of  global  nodes. 

NE  *  number  of  elements. 

n  ■  number  of  unknowns  in  the  structure. 

m  -  number  of  retained  modes  (if  no  reduction 
then  m*n). 

PathS  -  string  variable  indicating  the  choosen 
method  of  solution. 

{MSI )  and  (MS2)  provides  storage  for  the  mode 
shape  vectors.  The  1  and  2  refer  to  the 
and  Improved  Iterative  values. 

[li]  provides  storage  for  the  mass  matrix, 
its  dimensions  are  mxm. 

[K]  provides  temporary  storage  of  the 

stiffness  matrix  or  its  Inverse,  depending 
on  the  stage  of  the  program. 

[D]  provides  storage  for  the  result  of 
[Klinverse  *  [Ml.  it  is  used  to  iterate 
toward  the  correct  mode  shape. 

(TJ  if  the  structure  has  been  reduced,  this 
provides  storage  transformation  matrix 
to  convert  reduced  mode  shapes  to  full 
ones. 

IP  I  ]  and  {P2l  provide  storage  for  matrices 
used  in  constructing  the  above  trans- 
mation  matrix,  IT].  They  are  dimensioned 
(n-mtom. 

[Si  provides  storage  for  the  sweeping  matrix. 
This  matrix  is  used  to  remove  last  mode 
shape. 


Available  Sub-Programs: 

Display  .Matrix 
Store  .Matrix 
Retrieve  .Matrix 
Mat  .times  .Mat 
MatT r  arts  .times  .Mat 
Mat  .Plus  .Mat 
Invert  .Mat 


rnmmm  net  meet  frm 

IWUftflQ  RNUIH 


frvwfDJ 


u*. 

CDPlntP 

' 

OpFffpxwpli  *olv* 
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'  4 - 4 

'  I  OynFEP. eigen  solver  ! 

'  4 - 4 


COttON  GN,NE ,D0F ,n  ,n  ,PN$ ,Path* :  GOTO  Start 
'  Subroutines  Below 


Renove. Last .Mode: 

FOR  i=Al  TO  n:  Tll<Al,i>=0:  FOR  j=A!  TO  n:  T2»(i,j)=0:  NEXT  j:  NEXT  i:  T3KA1,A1)30  'init  temp  storage 

CALL  MatTrans. t ines .Nat (A1  ,n,n ,HSlt< )  ,NI< )  ,T1 8< ) )  '  [Til  =  {NSDtrantM 

CALL  Mat. tines.Mat(n,n,Al ,NS1I<) ,1180,1210)  '  IT2I  =  {MSntMSJtrantM 

CALL  Nat .  t  ines  .Mat  (A1  ,A1  ,n  ,T1  •<)  ,HS1 1( )  ,T38< ) )  7  T3  3  CMSUtranCHTCHSl } 

aa#=l/T3l(Al,Al):  CALL  Nat.plus.Mat<n,n,al,S»0,-aai,T2IO) 

'CALL  Retrieve .Matrix(n,n,T2*<), PNI4*. O’, A3)  '  load  original  EDI  into  T2# 

FOR  i=Al  TOn:  FOR  j=Al  TO  n:  T2i< i , j >=D#< i , j) :  0l(i,j)=0:  NEXT  j,i  '  init  ID) 

CALL  Nat .  t  ines  .Mat  (n  ,n  ,n  ,T2i< ) ,  SIO  ,DIO )  '  newCD]  3  originaHDHSIlatest 
RETURN 

Create. Eigen. Files: 

RL=n*8:  F<=PN$+’.S’:  OPEN  F*  AS  «1  LEN=RL:  F1EL0I1.RL  AS  BBi  '  the  shape  tile 
a$=MKW<0):  FOR  i=l  TO  n:  b$=M+a$:  NEXT  i  '  load  shape  tile  with  zeros 
FOR  i=l  TO  n:  LSET  BB4=bS:  PUTll.i:  NEXT  i 

RL-8:  F$=fW+’.  eigen’:  OPBt  F$  AS  *2  LB*=RL:  F1ELD»2,RL  AS  CM  '  tile  tor  eigenvalues 
RETURN 


Start:  Accuracy3. 01/100:  al*l:  a2=2:  a3=  3:  a*=l:  G0SU8  Create. Eigen. Files 
CALL  TEXTFONT(l):  CALL  TEXTSIZE(9):  IF  nOn  THEN  Flag*3’*’  ELSE  FlagV=** 

DIN  DNn ,n) ,MSl*(n ,A1 ) ,NS2«<n ,A1 ) ,T1 KA1 ,n) ,T3«<A1 ,A1 ) 

DIM  K*(n,n+1):1F  Flag!3’*’  THEN  CALL  Retrieve .NatrixIn.n.KIO.PW+’.Ke’.aS)  ELSE  CALL  Retrieve .Matrix(n,m4i, KUO, PNS4-.K 
&F’,a3) 

'liilHIlHiHHIiiHHIHHlIIIHHHH  'debug 

PRINT  USING  *n=  (ft  »=  ««  Flag*3  >!<’in,n,Tlagft 

IF  Flag*3’*’  THEN  CALL  Display  .Matrix<n,n,KftO,’Reduced  Stittness’)  ELSE  CALL  Display .Matr i x(i*,n4 1  ,KN< ) , ’St i ttness’ ) 
CALL  lnvert.Matrix(n,KftO) 

DIN  N#(n,n):  CALL  Retrieve .Natrix<n,n,N«0,PN«e,.M’eFlagf  ,a3) 

'  IHliiHHHIHIHHHilillHHIlHHIl  'debug 

CALL  Display.Natr ix<nTn-»l ,KttO .'Inverted  Stittness*) 

CALL  D i spl ay .Na tr i x (n ,n ,N8< ) , 'Mass’ ) 

'limHHHHIIIHHHIHHHIHHHIH 

CALL  Nat . t ines .Nat (n,n,n ,KI( ) ,N»< ) ,DI< ) )  '  ID3  3  IKIinvtNI 
CALL  Store .Matrix(n,n, DIO, PNO+’.D’jaS)  '  tenporary  tile 
ERASE  K8  '  clear  sane  nenory 

DIN  SI(n,n),T2l(n,n):  FOR  i^Al  TOn:  S«(i ,i)=Al:  NEXT  i  'init  (SI  as  III 

IF  nOn  THEN  DIN  Pl»<n,n),P2ft<n,m):  CALL  Retrieve.Natrix<n-ft,n,Pl»0,PN$*’.Pr,a3):  CALL  Retrieve .Natrix(n-n,n,P2»0,PNft 
♦*.P2’,a3)'load  PI  It  P2 

FOR  eigen=Al  TO  n  '  begin  solution 
FOR  Nil  TO  n:  MS1»( i  ,A1)=A1 :  NEXT  i 

i=2:  j=eigen:  Sign=-Al  '  create  a  tirst  guess,  should  have  one  less  sign  change  as  eigen 
WHILE  j)=i:  FOR  k=i  TO  n:  NSIKi ,Al)=Sign:  NEXT  k:  Sign=-Sign:  i=i*Al :  WEND 

Change=Al 
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WHILE  Chang*>Accuracy 

CALL  Hat. tints. Mat (i»,Al ,n,D«( ) , MSI  10 ,MS2I( ) ) 

Freq2*=l/HS2»<A1,A1):  Changed 
FOR  i=Al  TO  n 

MS2K  i  ,A1)=MS2»<  i  ,Al)*Freq2# 

NunFABS<  (MS1K  i  ,A1  )-MS2V<  i  ,A1 ) )/HS2l(  i  ,A1 ) ) :  IF  Nun>Change  THEN  Change=Nun 
MSI »( i ,A1 )=MS2I< i ,A1):  MS2»< i ,A1 )=0 
NEXT  i 

'  HiHHHHHHIHHHHHHIHHHiHl  '  debug 

IF  2=0  THEN  CLS 

CALL  M0WET0<2,50):z=z*l:  PRINT  USING  ,Try*ll,;2 
PRINT  Change 

'CALL  Display .Natrix<n,al,NSll<), ’Trial  Vector1) 

'  HHHIHHHHIHHHHHHHIHHHH 

WEND: z=0 

IF  n=m  THEN  FOR  i=l  TO  n:  MS2l(i,Al)=MSll(i,Al):  NEXT  i:  GOTO  Skip  'transtom  it  structure  not  reduced 
FOR  i=l  TO  n-i»:  FOR  j=l  TO  n:  Tl(i,j)=0:  NEXT  j,i  '  store  CPU  in  IT] 

CALL  Mat . plus. Mat(n-n,n,al, TOO ,Freq2l,P2l<))  '  tind  m  then  below  create  natrix  with  III  over  -IT) 
FOR  i=l  TO  n-n:FOR  j=l  TO  i»:TI<i+n,j)=-TKi,j):NEXT  j.isFOR  i=l  TO  i»:FOR  j=l  TO  n:T*<i,j)=-<i=j):NEXT  j,i 
FOR  i=n  TO  n:  MS2K i ,A1  )=0 :  NEXT  i  'Finish  initialling  MS2I 
CALL  nat . t ines .Mat <  n ,A1 ,n ,TI( ) ,MS1 8< ) ,MS2I( ) )  '  Full  eigenvector  stored  in  MS2# 

Skip: 

'  HHHHHHIHHWHHIHIIHHHHH  '  debug 

PRINT  USING  ‘Mode  Shape  ««  Freq  Ml.ll . ;e i gen ,SQR<Freq2#> 

CALL  Display .Natrix<n,Al ,MS2»<),‘ Mode  Shape1) 

'HHiHHiHHHHHHiHHHIHHiHH 

• 

'Prepare  tor  next  node  shape 
Lt=(eigen-Al)*8:  Rt=<n-eigen)»8 
FOR  i=l  TO  n:  GET41 , i :  at=BB$ 
a$=LEFTf <  a< ,Lt )  +MKM<MS2*<  i ,  1 ) ) +RI GWTV  a$ ,  Rt ) 

LSET  BB*=a*:  PUTil , i  '  store  elenent  in  eigen  vector  natrix 
NEXT  i 

LSET  CC$=MKDf<Freq2l):  PUTI2, eigen  '  store  square  ot  eigenvalue 
IF  eigentn  THEN  GOSUB  Renove. Last  .Mode 
NEXT  eigen 

CLOSE:  KILL  PtW.D1:  KILL  PN$4*.k»*:  KILL  PW+'.Ne1:  KILL  PNH'.Pi1:  KILL  PN$*1.P21  '  distroy  tenporary  tiles 

'HiiHIHHHmHHHHIHHHIHtlHI  'debug 

CALL  Retrieve .Matrix(n,n,S«<), Pm*1. S'.al) 

GALL  Display .Matrix<n,n,S«), ‘Mode  Shapes1) 

CALL  Retrieve .Matrix(n,al ,EI<) ,PN$+’ .eigen* ,al) 

CALL  Display .Matr i x<ntnl ,Ei() ,'Eigen  values1) 

'HiKHifilliKHIliHlillKilKHillil 

CHAIN  ’Basic  Disk  lsDynFEP. uncouple/solve1 
END 


Sub-Prograns  Below 
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I  tutlc  fa 
(Ft 


DynFEP.uncoupfe/solve 

_ Flow  Diagram _ 


Notes: 


Flnd/Stera  gmaralizad  imm  matrix 

I _ WiWjg _ 


Program  Variables: 


FMOhnvana 


load  Initial 
CooOltlona  [Ual 


Apply  Boundary  Conditio 


Find  panaralltad  ca 

BJol .  lailUel 


land  boundary  conditio 

_  (BC) 


Structuri 

reduced? 


load  maaa  matrix 
W 


laad  atganvoluao 

to 


Find  Nawmark  conotanlo 
A  QtC,  A3.  Ad.  A7 


Find  dynamic  forcaa  from 
forcaa  an  nodaa,  forcaa  on 
alamonto.  Inorttal  forcaa 
from  aaaanttal  BC 


Find  gonarotlnd  force  matrix 
_ IFd)  ■  IStKffd)  ♦  IFoll 


9g|yg  goOF  I^UItlOflt 

Ull  a  (rot mi  ♦  AO*UOI  ♦  A2-UOI  a  A3*U0D  /  (AO  *  Cl) 


Find  currant  ganarallzod  oceorotlon  t  velocity 
wn  a  A1*((UI)  -  (UOI)  -  AS-WO)  -  A3»(U0M 
tun  a  (uoi  •  at-  oai>*njo-)  ♦  oai*nin)n>oitaT 


Swap  nil)  and  (UO> 
for  next  tuna  otap 


FMSaioro  real  caord.  valuao 

wi,  ur.uii  a  i3iui,  ur.  un 


GN  *  number  of  global  nodes. 

NE  -  number  of  elements. 

n  *  number  of  unknowns  in  the  structure. 

m  -  number  of  retained  modes  (If  no  reduction 
then  m*n). 

Path$  -  string  variable  indicating  the  choosen 
method  of  solution. 

[S]  provides  storage  for  the  mode  shape 
matrix. 

[MI  provides  storage  for  the  real  or  generalized 
mass  matrix,  depending  on  the  stage  in  the 
program.  Its  dimensions  are  mxm. 

(UOI  and  [U 1  ]  provide  storage  for  the  current 
and  last  generalized  displacement, 
velocity,  and  acceleration  vectors. 

Treated  mathematically  as  3  column 
matrices,  they  are  stored  as  a  nx3. 

[U1  provides  storage  for  the  current  real 
displacements.  Its  dimensions  is  nxl. 

[E]  provides  storage  for  the  eigenvalue 
matrix.  Mathematically  It’s  a  square 
diagonal  matrix,  it's  stored  as  a  mx  1 . 

(FdJ  provides  storage  for  generalized  dynamic 
forces.  It  is  dimensioned  as  mx  1 . 

(FsJ  provides  storage  for  static  forces.  It  is 
dimensioned  as  nx  1 . 


Available  Sub-Programs: 

Display  .Matrix 
Store  .Matrix 
Retrieve  .Matrix 
Mat  .times  .Mat 
MatTrans.  times  .Mat 
Mat. Plus  .Mat 
lnvert.Mat 


Ftfriatod? 


|  Htcrumrt  Time 

extern  te 

DynTir  .■  l  i » \  ■ 
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'  4 - 4 

'  I  DynFEP. uncouple/ solve  I 

'  4 - 4 


CQttl ON  GN ,NE, DOF, n,m,PN$, Path!:  GOTO  Start 
'  Subroutines  Below 


ReadHistoryFile: 

OPEN  History!  AS  14  LB*=ld:  FIELD  14,8  AS  2!<1),8  AS  2!<2) 

GET»4,1:  Max=CVS<2!<D):  Min=2:  i=lNT<<Max-Min)/24l) 

6ET»4,i:  T1=CVS<2!<1)) 

WHILE  Nax)Min4l  AND  TlOT+phaie 
IF  Tl<T4phaze  THEN  Nin=i:  i=Nin*!NT«Nax-i)/2) 

IF  TDTephaze  TH04Max=i:  i=Max-INT«i-Nin)/2) 

GETi4, i :  T1=CVS<2!<1)) 

WEND:  IF  T4phaze=Tl  THEN  4=CVS<2!<2)):  60T0  Found 
6ET84 ,liin:  T1=CVS<2!<1)):  +  1=CVS<2S<2)> :  IF  T+phaze=Tl  THEN  4=41:  GOTO  4ound 
6ETI4,Max:  T2=CVS<2!<1)):  42=CVS<2!<2)):  IF  T*phaze=Tl  THEN  4=42:  GOTO  4ound 
4=<T4phaze-Tl)«(42-41)/(T2-Tl)441  '  interpolate 

4ound:  CL0SEI4:  RETURN 

AssenbleForceNat:  1ast=0 
FOR  i=l  TO  GN 

FOR  j=0  TO  DOF-1:  index=< i-l)*DOF4j4l 
IF  8CK index,  1)=1  THB4  GOSUB  NodeDynForces 
NEXT  j , i :  RETURN 

NodeDynForces:  k=345*j:  DloacH):  A=l4j»4:  IF  lastOi  THEN  6ET9] , i :  last=i  'set  index's  and  read  node  dynamic  4orces 
IF  MIM(Flgl!,A42,2)=,ir  THEN  anp=CVS<N!<k4l)):angle=CVS<N!<k42)):phase=CVS(N!<k43)):Dload=amp»SIN<t*angle4phase)'Harm 
Load 

IF  MID9<FlglS THEN  History!=E!<k44):G0SUB  ReadHistoryFile:Dload=4  '  non-harmonic  load 
Fdl(index,l)=Fdl<index,l)4Dload  '  add  in  Force;  positive  to  right,  upward,  clockwise 
RETURN 

Inert ialForces: 

FOR  i=l  TO  2:  6ETil ,EBCZ< i ,1):  j=EBC«i,2):  AF3+<j-l)*4:  k=3+(j-l)*5:  Displ=0 

IF  NID4<FlgU,A*2,2)=,ir  THEN  a»p=CVS(N$(k+l)):angle=CVS(N«(k42)):phase=CVS<N$(k*3)):Displ=DispHamp»SlN(t*angle4phas 

e) 

IF  NIIMIFlgli^^^'lO'  THEN  History!=E!(k44):60SUB  ReadHistoryFile:Displ=DispH4  '  non-harmonic  Displacement 
FOR  k=l  TO  GN:  index=(k-l)#DOFM:  Fd#<index,l)=-Displ:  NEXT  k 
NEXT  i:  k=0 

FOR  i=l  TO  GN*NE  '  apply  essential  BC 

IF  BCI(i,l)=l  TH04  k=k4l:  IF  iOk  TH04  SUAP  Fdll<k,l),Fd»<i ,1) 

NEXT  i 

CALL  Mat . t imes .Nat (n ,n 1 , n ,NI< ) ,Fdl< ) ,T2I< ) ) 

RETURN 

El etaen  tMatr i xAssenbl er : 

FOR  ELEMENT2 1  TO  NE 
60SUB  BuildEI emen tMatr ices 
IR=<N1X-1 )»D0F : IC=(N2X-1 ) eDOF 
' —  Asseii.Dyn.  Elen.  Forces: 

FOR  i>l  TO  DOF:  Fd*<]R4i,l)=Fd«<IR4i,l)44«(i):  Fdl(IC4i,l)=Fd«IC4i,l)44l(i4D0F):  NEXT  i 
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BuildElenentHatrices: 

GETI2 , ELENENT :N  15NCVI  <LtS)  :N2K=CV  I  <  R  tf )  '  get  left  and  right  global  node  Is 
GETf 1 ,N17. :X1=CVS(N$< 1 ) ) :Y1=CVS(N$< 2) )  '  get  left  side  coord's 
6ETI1 ,N2X:X2sCVS(NJ(l)):Y2sCVS(N$(2))  '  get  right  side  coord's 
L=SQR«Y1-Y2)*2+(X1-X2)*2)  '  find  elenent  length 
' —  Elen. Dyn. Forces: 

DDloa<NO:TDload=Q 

IF  MIW<Fl92»,2t2)=Ml'  THEN  anp=CVS(E*<10)):angle=CVS<Ei<ll>):pha2e=CVS<Ei<12)>:DD1oad=anp*SIN<t*angle+phase)'Harn 
Dyn 

IF  NIW(F1g2$,2,2)s‘10‘  THEN  History$=EI(13):60SUB  ReadHistoryFile:DDload=f  '  non-harnonic  dyn  load 
IF  NID$<Flg2$,5,2)=,ll‘  THEN  an^S(E«ld)):angle=(^S<E«(17)):pha2e-(MS(E$(18)):TDload=anp>SIN(t>angleiphase)'Harn 
TanDyn 

IF  MID*<F1g2$,5,2)=,,10“  THEN  History4=E*(19):60SUB  ReadHistoryFile:TDIoaiNf  'non-harnonic  tan  dyn  load 
f  l< 1 )=L*TD1 oad/2 :f l( 4)=f l< 1 )  '  positive  to  the  right 

fl(2)sL*DDload/2:fl<5)=ft(2)  '  positive  upward 
fl(3)=-L‘2*DDload/12:fl<4)=-f#<3)  '  positive  clockwise 

IF  Xl-X2=0  THEN  angle=SGN<Yl-Y2)*Pil/2  ELSE  angle=2»Pi«-ATN«Yl-Y2)/(Xl-X2)) 

IF  ang1e>2»Pi#+.003  OR  angle<2»Pil-.003  THEN  TransfornDynForce 
RETURN 

TransfornDynForce:  'Subroutine  to  transforn  Stiffness,  Nass,  and  Force  elenent  matrices 
' —  BuildTransfomationMat:  '  build  IT] 

FOR  i*l  TO  6:F0R  j=l  TO  d:TI(i ,j)=0:NEXT  j,i  '  initialue 

T«( 1 , 1 )=C0S<  angle) :T«<4 ,4)=T*< 1,1) :T*<  2 , 2)=T«< 1,1) :T»(  5,5)=T*< 1,1) 

TI<1 ,2)=-SIN<  angle)  :TK4,5)=Ti(  1,2)  :TK2,l)a-TI<l,2):TI<5,4)=TI(2,l) 

TI(3,3)=1:T*(6,4)=1 

FOR  i-1  TO  <S:CKi,l)=0:F0R  k=l  TO  4:CI(i,l)*C»(i,l)*T«<k,i)*f»(k):NEXT  K,i:FQR  i=l  TO  d:fl(i)=C»(i,l):NEXT  'tran[T)*{f) 
RETURN 

Get .del taT.and.T ine .Steps: 

Del taT=2/S0R(E»(n,  1 ) ) :  T4=STR«DeltaT):  TS#=STRS< INT<SQR<ER<1 , 1  ))/Del taT)*l )  '  nax  del taT  and  nin  Icycles 
UINDOU  3,,<250,22)-<505,132),-4 

CALL  TEXTSIZE(IO):  CALL  MOVETO(5,2d):  PRINT  ‘Enter  tine  step  (nax.  shown)*:  CALL  TEXTSIZEC12) 

EDIT  FIELD  l,T$,(5,30)-(250,45) 

CALL  TEXTSIZE(IO):  CALL  MOUETO<5,<1):  PRINT  ‘How  nany  tine  steps?’:  CALL  TEXTSIZEU2) 

EDIT  FIELD  2,TS$,(5,d5)-<250,80>:  EDIT  FIELD  1 

BI/nON  1 ,1  ,‘OK*  ,<200 ,84)-<250,102) 

i=l 

loop: 

<tsDIAL06<0) 

IF  #1  OR  iNd  TH04  done  'got  OK  button  or  RETURN 
IF  <N2  THEN  i=01AL0G<2):  EDIT  FIELD  i  'got  field  selection 
IF  <N7  THEN  i=(i  NOD  2)+l:  EDIT  FIELD  i  'got  TAB  key 
6OT0  loop 

done:  CALL  TEXTSIZE(IO):  DeltaT=VAL(EDlT$<D):  NunSteps=WL<EDIT*<2)):  UINDOU  CLOSE  3 
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FIRD«1,12  AS  FlglS,  4  AS  N$<1),  4  AS  N$<2),  4  AS  N$<3),  4  AS  NS<4>,  4  AS  Nf(5),  4  AS  8  AS  N$<7),  4  AS  N$<8),  4  A 
S  NS<9),  4  AS  N$<10),  4  AS  NS<11),  8  AS  N$<12),  4  AS  N$<13),  4  AS  N$<14),  4  AS  Ni<15),  4  AS  N»<16),  8  AS  N$<1?) 
F<=PN$e*.elenents*:OP0f  Ft  AS  12  LEN=94 

F1ELD*2,6  AS  F1g2$,2  AS  L»,2  AS  Rt*,4  AS  E*<1),4  AS  Et<2),4  AS  Et(3),4  AS  Et<4),4  AS  E*<5),4  AS  Et<4),4  AS  E*<7),4  AS  E 
1(8), 4  AS  E*<9),4  AS  E$<10),4  AS  E*<11),4  AS  E«12),8  AS  Et<13>,4  AS  E$<14),4  AS  E$<15>,4  AS  E*<14),4  AS  E$<17),4  AS  Et( 
18)  ,8  AS  E«19) 

F4=PN44\displ*:  OPEN  Ft  AS  43  LEN*24 
FIELD#3,8  AS  Ut,  8  AS  Vt,  8  AS  Act 

IF  n<)n  THEN  Flag$='«*  ELSE  Flag$="  '  flag  reduced  structure 
DIN  KI(n,n4l),Fs8<n,l):  CALL  Retrieve  .Natrix<n,n4l,KIO,PN$4*.K4F,,n4> 

FOR  i=l  TO  n:  FsA< i ,1)=K»< i ,n+l) :  NEXT  i:  ERASE  Kl  '  load  static  force  matrix 

'Find/Store  General i zed  Nass  Natrix 

DIN  NI(n,n),Ndia«ii,nl),Tll<n,n)i  CALL  Retrieve .Natrixtn.n.NHO.PNt+’.N'.nt) 

DIN  Sl(n,n):  CALL  Retrieve.Natrix<n,n,SI(),PN$4*.S*,n4)  '  load  node  shapes 
CALL  NatTrans. t ines .Mat (n ,n , n , Sl< ) ,N»< ) ,T1*< ) ) 

FOR  i-\  TO  n:  FOR  j=l  TO  n:  Nl< i ,j)-0:  NEXT  j,i  'init  HI 

CALL  Nat.tines.Nat(n,n,n,Tll<)fSIO,N*<)):  FOR  i~l  TO  n:  Hdia«(i,l)=NI<i,i):  NEXT  i  '  store  dia.  in  Ndial 

'  HHHHiHHHHHHHHHHHHiH  debug 

CALL  Display .Hatrix<n,n,N»<), ’Generalized  Nass  Natrix*) 

CALL  Display  .Matrix(n, n, SIO , ’Node  Shapes*) 

'iHlHHHHilllHlilHHilllHIHII 

CALL  Retri eve. Natr ix(n,n,N*<),PN$4*.N*,n4)  '  INI  needed  to  find  inertial  forces 
ERASE  TIP  '  clear  sene  nenory 

DIN  El(n,l):  CALL  Retrieve.Natrixta^^EK^PNt+'.eigen'tnt)  '  load  eigenvalues 

DIN  8CI< GN*D0F ,  1 ) : p=GN*D0F :  CALL  retri  eve. Natr  ix<p,nl,BC»0,PN$4*.BC*,n4)  '  load  boundary  condition  index 

'  HHHHIHHIHHHHIHHHIHHH  debug 

CALL  Display .Natrix<GN»DOF,nl,BCI<), ’Boundary  Condition  Index*) 

'iHHHHiHHiHHHlHIHHIHHH 

DIN  U0l(n,3) ,U1N<GN«NE ,3) :  GOTO  Skip  'trouble  with  initial  conditions  file,  can't  resolve 
IF  Flagts***  THEN  Skip  '  initial  conditions  nust=0  if  structure  is  reduced 
CALL  Invert.Natrix<n,S«)) 

'  »*»«**M*******m***m»***m**M  debug 

CALL  Display .Natr ix(n ,n ,SMO , 'Inverted  Node  Shapes') 

CALL  Retrieve .Matrix<GN«00Fln3(Uli<> ,PNt4* .initial*, n4):  k=0 
FOR  i=l  TO  GN*D0F  'apply  boundary  conditions  to  intial  conditions 
IF  BCI(i,l)=l  THEN  Mel:  IF  iOk  THEN  FOR  j-1  TO  3!  SWAP  01«(k,j),01*(i,j):  NEXT  j 
NEXT  i 

CALL  Nat. tii»es.Nat(n,nl,n,SK),Ul»O,U0»())  '  find  generalized  initial  conditions 
CALL  Retrieve.Natrix(n,n,SI(),PN$4*.S*,n4)  '  reload  node  shapes 
Skip:  'CALL  Store .Natr ix<n ,n3,UlR< ) ,PNSes .di spl * ,n4)  'store  initial  conditions  in  displacenent  file 

Find. Essential .6C:  i=0:  jsi 

UH1LE  j<=2:  i-i *1 :  index=<i-l)»D0F4j  '  find  out  where  unifom  base  novenent  stored 
IF  8C«< index, l)-0  THEN  EBCX(j,l)=i:  EBCX(j,2)=j:  j=j+1  :  i=0 
WEND 

'6et  or  calculate  constants 
delt)Fl/2:alpha?l/4:  GOSUB  6et.deltaT.and.Tine. Steps 

Atel/<alpha»DeltaT*2):A2»l/<alpha*0eltaT):A3*l/<2»alpha)-l  '  calculate  constants 
Ad>Del taTf < 1 -del ta)  :A7»del tafDel taT 
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FOR  Countcr-l  TO  NumSteps  '  begin  solution  loop 

DynForceNat:  FOR  i=l  TO  QN«NE:  Fd«<i,l)=0:  NEXT  is  FOR  i=l  TO  ns  T2t<i,l>=0:  NEXT  i  'init  Fdl  &  T2* 

60SUB  InertialForcess  GOSUB  AssembleForceMat:  60SUB  ElementNatrixAssembler 
FOR  i=l  TO  GN#F£  '  apply  essential  BC 

IF  BC»<i,l)=l  THEN  k=K+ls  IF  i<)k  THBF  SNAP  Fdl<K , 1 ) ,Fd8< i , 1 ) 

NEXT  i 

CALL  Hat. P1us.Nat<n,nl,onel,FdK),onel,T28<))  'add  node/elenent  Forces  and  inertial  Forces 
CALL  Mat.Plus.Mat(n,nl,onet,Fdll<),onel,Fsl<))  'add  dynamic  and  static  Forces 
FOR  i=l  TO  ns  T2«(i,l)=0:  NEXT  i  'init  T2» 

CALL  MatTrans.times.Mat(n,nl,n,S80,Fdf0,T2#0)  '  Find  generalized  dyn.  Force  matrix 
Solve: 

FOR  i=l  TO  m 

Ul«( i ,l)=(Fd»<i ,1)/Mdia»< i , l)+AO*UOi< i ,1)*A2*U0I( i ,2)+A3*U0#< i ,3) )/<A0+E#< i ,1)) 

NEXT  i 

FOR  i=l  TO  m  '  Find  V  and  A  vectors  and  store  displacements  in  U 
Ul»<  i  ,3)=A0*(U1II(  i  ,1)-U0»<  i ,  1 ) ) -A2*U0»<  i  ,2)-A3*U0»(  i  ,3) 

U1K  i  ,2)=U0»<  i  ,2)*A6*U0«<  i  ,3)+A7*Ul«(  i  ,3) 

FOR  j=l  TO  3:  U0«< i ,j)=Ul«<i ,j)s  Ul*<i,j)=0:  NEXT  j  '  IU0I  =  CU11  For  next  time  step,  init  IU13 
NEXT  i 

Find. Store. real. displacements: 

CALL  Mat. times.Mat(n,n3,m,SI(),  1)010,11110)  '  (Ull  =  ISHUOI  convert  From  general ized  coordinates 

' *******»****}**************»****  debug  only 

CALL  Display .Natrix<n,n3,Ul»<), ‘Displacement,  Velocity,  Acceleration*) 

'HKHHHHHHHHHHHIIHH 

FOR  i=l  TO  n 

LSET  W=HKSttUll<i,l)):LSET  V*=MKS*<Ul*<i ,2)):LSET  Ac<=HKSl(Ull(i ,3)):j=i*Counter*n:PlfT»3,j  '  save  to  disk 
NEXT  i 

T=T+De1taT  'next  time  step 
NEXT  Counter 

CLOSE:  CHAIN  'Basic  Disk  l:OynFEP.menu* 

BFD 
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Program  Variables: 

6N  -  number  of  global  nodes. 

NE  =  number  of  elements. 

n  -  number  of  unknowns  In  the  structure. 

m  =  number  of  retained  modes  (if  no  reduction 
then  m-n). 

Path$  -  string  variable  indicating  the  choosen 
method  of  solution. 

[K]  provides  storage  for  the  global  stiffness 
matrix,  and  the  static  force  matrix. 

The  stiffness  matrix  is  a  square  matrix 
with  dimensions  equal  to  6N  times  3. 

The  static  force  matrix  is  stored  with 
the  stlffhess  in  an  addition  column. 

[Ml  provides  storage  for  the  global  mass 
matrix.  It  is  a  square  matrix  with 
dimensions  equal  to  GN  times  3. 

[UOJ  and  [U 1 J  provide  storage  for  the  current 
and  last  displacement,  velocity,  and 
acceleration  vectors.  Treated  a?  3  column 
matrices,  they  are  stored  as  a  nx3. 

fFd]  provides  storage  for  generalized  dynamic 
forces.  It  is  dimensioned  as  mx  1 . 


Available  Sub-Programs: 

Display  .Matrix 
Store  .Matrix 
Retrieve  .Matrix 
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'  ♦ - ♦ 

'  I  OynFEP  I 

'  ♦ - ♦ 


COtlQN  GN ,NE, DOF, n,n,PN), Path):  GOTO  Start 
'  Subroutines  Below 


ReadHistor/Fi le:  BUTTON  10,2 

OPEN  History)  AS  14  LEN=14:  FIELD  14,8  AS  2X1) ,8  AS  2X2) 

GETI4.1:  Nax=CVS(2Xl»:  Min=2:  i=lNT«Max-Nin)/2+l) 

6ET#4,  i ;  T1=CVS(2XD) 

UHILE  Nax)Min+l  AND  TIOT+phaze 
IF  THT+phaze  THEN  Nin=i:  i=Min*lMT<<Max-i)/2) 

IF  Tl)T+phaze  THEN  Max=i:  i=Max-lNT<<i-Nin)/2) 

6ET*4,i:  T1=CVS(2XD) 

WEND:  IF  T+phaze=Tl  THEN  f=CVS(2X2)):  GOTO  found 
GETI4,Nin:  T1=CUS(2XD):  fl=CVS(2X2)):  IF  T+phaze=Tl  THEN  f=fl:  GOTO  found 
GETI4,Nax:  T2=CVS(2XD):  f2=CVS(2X2)):  IF  T+phaze=Tl  THEN  f=f2s  GOTO  found 
f=(T+phaze-Tl)»(f2-fl)/(T2-Tl)+fl  '  interpolate 
found:  CL0SEI4:  BUTTON  10,1:  RETURN 

Guass:  BUTTON  14,2 

FOR  i=l  TO  n:N#=KXi,i):FOR  j=I  TO  N+l:KXi,j)=K»<i,j)/N»:NEXT  j 

FOR  k=l  TO  n:lF  kOi  THEN  M*=K*(k,i):F0R  j=i  TO  n+l:K#<k,j)=KI(k,j)-KKi ,j)*M»:NEXT  j 

NEXT  k,i:  BUTTON  14,1:  RETURN 

AsswibleForceNat:  BUTTON  7,2 
FOR  i=l  TO  GN;  GET*l,i  '  read  specified  nodal  loads 
FOR  j=0  TO  DOF-1 
index=< i-l)«DOF+j*l :k=3*5*j 
' —  NodeDynForces: 

DloaiH):  Af3*j*4 

IF  MIDXFlgl),A,2)=Mr  THEN  ar»p=CWS(N)(k*l)):angle=CVS(NXk42)):phase=CVS<NXk43)):Dload=anip*SIN(t*angle*phase)'ha 
rn 

IF  NIDXFlgl),A,2)=,10*  THEN  HistoryXEXk+4):G0SUB  ReadHistoryFile:Dload=f  '  non-harmonic  load 
Fdt<index)sFdl<index)*Dload  '  add  in  force;  positive  to  right,  upward,  clockwise 

Ul#( index ,1)=0 :k=l+4*j : IF  MIM<F1  glS ,k ,  1 THEN  Ultt< index ,1  >-l  '  Flag  essential  B.C.,  used  later 
NEXT  j,i:  BUTTON  7,1:  RETURN 

Essential .B.C:  BUTTON  13,2 

Disp1=0:  Node=lNT((i+2)/DOF):  j=< i+2)  NOD  DOF:  k=3+j»5:  A=H4*j:IF  j=0  THEN  GET II  ,Node 
IF  NJDXF!gl$,A*l,l)=M*  THEN  Oispl=CVS<NXk))  'static  displacement 

IF  MlDXFlgl),A42,2)=,ll*  THEN  anp=CVS<N)(k*l)):angle=CVS<NXk*2)):phase=CVS(NXk*3)):Displ=Displ4amp*SlN(t*angle«phase 

) 

IF  NlDXFlgl),A42,2)=*10*  THEN  History)=EXk44):G0SUB  ReadHistoryFile:Displ=Displ4f  '  non-harmonic  Displacement 
BUnON  13,1:  RETURN 

ElenentNatrixAssembler:  BUTTON  8,2 
FOR  ELEMENT5 1  TO  NE 
GOSUB  BuildElementNatrices 
I R*(N1Z-1 ) *D0F : I C=(N2X-1 ) »D0F 
' —  Assem.Dyn.Elem. Forces 

FOR  i*l  TO  DOF:  Fd»< IR4 i )=Fd*< IR+ i )4f #( i ) j  Fdf< 1 €4 i )*Fd»< I C* i ) *f •< i ♦ DOF) :  NEXT  i 
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NEXT  elenent:  BUTTON  8,1:  RETURN 
Bui IdElenentNa trices:  BUTTON  9,2 

GETI2, ELEMENT :NlX=CVI(Lt*):NZ<=CVI(Rt*)  '  get  left  and  right  global  node  »'s 
GET81 ,N1X :X1=CVS(N$( 1 ) ) :Y1=CVS(N$( 2) )  '  get  left  side  coord's 
GETI1  ,N2Z :X2=CVS(N$( 1 ) ) :Y2=CVS(N$(2) )  '  get  right  side  coord's 
L=SQR((Y1-Y2)‘24(X1-X2)‘2)  '  find  elenent  length 
' —  Elen. Dyn. Forces: 

DDIoacM  :TDload=0 

IF  MlDi<F1g2$,2,2)='H'  THEN  anp=CUS<E$<10>):ang1e=CyS<E*(ll)):phaze=CVS(E$<12)):DDload=anp*SIN<t*angle4phase)'harn 

IF  MID$(Flg2$,2,2)='10'  THEN  History*=E$< 13): GOSUB  ReadHistoryFile:DDload=f  '  non-harnonic  dyn  load 

IF  MIW(Flg»,5,2)=*ll*  THEN  anp=CVS(E$(14)):angle=CVS(E«17)):phaze=CVS<E$(18)):TDload=anp*SIN(t*angle4phase)'harn 

IF  MIDS<F^g2$(5,2)=M8,  THEN  History*=E$(19):G0SUB  ReadHistoryFile:TDload=f  'non-harnonic  tan  dyn  load 

ft(l)=L»TDload/2:fl!(4)=fll(l)  '  positive  to  the  right 

•flK2)=L*0Dloacl/2:-fll<5):='fi<2)  '  positive  upward 

f#(3)=-L‘2*D0load/12:fll(4)=-fil<3)  '  positive  clockwise 

IF  Xl-X2=0  THEN  angle=SGN(Yl-Y2)*Pi«/2  ELSE  ang1e=2»Pi«-ATN<(Yl-Y2)/<Xl-X2>) 

IF  angle>2*Pi 14.003  OR  angle(2»Pii-.003  THEN  GOSUB  TransfornDynForce 
BUTTON  9,1:  RETURN 

TransfornDynForce:  BUTTON  1 1 ,2: 'Subrout ine  to  transforn  Stiffness,  Mass,  and  Force  elenent  natrices 
GOSUB  BuildTransfornationMat 

FOR  i-1  TO  4:C«(i,l)=0:FOR  K=1  TO  6:C*<i  ,l)=C*<i ,  1  )+T*<l( ,  i  )*f  «(k>  .-NEXT  k,i:F0R  i=l  TO  6:f*(  i)=C*(  i  ,1)  :NEXT  'tran[T]*{f) 
BUTTON  11,1:  RETURN 

8uildTransfornationMat:  BUTTON  12,2:  '  build  IT! 

FOR  i=l  TO  4:F0R  j=l  TO  4:T»( i ,j)=0:NEXT  j,i  '  initialize 
IF  angle  MOD  Pill/2  THEN  T*(1,1)=C0S< angle)  ELSE  T*<1,1)=0 
T#<  4 ,4)=T*< 1 , 1 ) :T»<  2 , 2)=T»( 1 , 1 ) :T»<  5 , 5)=Ti< 1 , 1 ) 

IF  angle  MOO  Pi*  THEN  TI(l,2)=-SIN(angle)  ELSE  T»<1,2)=0 
T»( 4 , 5)=T»( 1 ,2) :T»< 2 , 1 )=-T»( 1 , 2) :TI( 5 ,4)=T»( 2 , 1 ) 

TI(3,3)=1:T»(6,6)=1 
BUnON  12,1:RETURN 

Get. del taT.and.Tine. Steps: 

T$=’Enter  tine  step  t.’:  TS$=*Now  nany  tine  steps?*  '  nax  del taT  and  nin  Rcydes 
UINOOU  3,, (250, 22) -(505, 132), -4:  CALL  TEXTFONT(l) 

CALL  TEXTSI2E(12):  CALL  M0VET0(5,24);  PRINT  'Enter  tine  step  (nax.  shown)':  CALL  TEXTS12E(12) 

EDIT  FIELD  l,T$,(5,30)-(250,45) 

CALL  TEXTSIZE(12):  CALL  M0VET0(5,4l):  PRINT  'How  nany  tine  steps?':  CALL  TEXTS1ZE(12) 

EDIT  FIELD  2,TS$,(5,65)-(250,80):  EDIT  FIELD  1 

BUTTON  1 , 1 , *0K* , ( 200 , 84) -<  250 , 1 02) 

i=l 

loop: 

d=DIALOG<0) 

IF  &=l  OR  d=4  THEN  done  'got  OK  button  or  RETURN 
IF  d=2  THEN  i=DIAL0G<2>:  EDIT  FIELD  i  'got  field  selection 
IF  d®7  THEN  i=(i  MOD  2)41:  EDIT  FIELD  i  'got  TAB  key 
GOTO  loop 

done:  CALL  TEXTF0NT(4):  CALL  TEXTSIZE<9):  DeltaT=VAL(EDIT«D):  NunSteps=MAL(EDIT«(2)):  WINDOU  CLOSE  3 
RETURN 

BigTextsCALL  TEXTFONT < 0 ) : CALL  TEXTSIZEt 12): RETURN  '  Chicago 
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LittleTextiCALL  TEXTFONT ( 1 )  :CALL  TEXTS12E(9) '.RETURN  '  Geneva 
NormalText:CALL  TEXTFONT <  1 ) : CALL  TEXTSIZE(IO)  jRETURN  '  Geneva 
Forma tedText: CALL  TEXTF0NT(4):CALL  TEXTS1ZE(  9)  .‘RETURN  '  Monaco 


Start:  CR1=CHR!<13):  Pi!M*ATN<l):  nl=l:  n$=3:  n4=4:  n=GN*OOF :  GOSUB  FormatedText 

'  create  status  windous 
F!='DynFEP.info*:OPEN  f!  AS  HI  LEtMO 

FlELDttl,  2  AS  XI!,  2  AS  YU,  2  AS  X2!,  2  AS  Y2!,  30  AS  Title!,  2  AS  Type! 

UINDOU  2,‘DynFEP  Input/Output  Window’, (14, dl)-(512,2d3),l 
UINDOU  1 ,’DynFEP  Status  Window’, <4, 41M424, 141), 1 
FOR  i=l  TO  1»:  GETttl ,  i 

xl=CVl(Xl!):  yl=CVl (Yl!) :  x2=CVl(X2!):  y2=CVl<Y2!):  A!=Title!:  K i nd=CVl (Type!) 

WHILE  RIGHT!(a!,l)=*  ■  :a!=LEFT!(a!,LEN<a!)-l) :W04D:BUTT(]N  i ,l,a!,<xl ,yl>-<x2,y2) .Kind 
NEXT  i : CLOSE! 1 

DIN  K*(n,n+l),Mt(n,n) ,Fd»(n) ,U0!l<n,3) ,Ull(n,3) ,N!<  17) ,E!< 19) 

F!=PN!+\ nodes'  :OPEN  F!  AS  #1  LEN=92 

FlELDil ,12  AS  Flgl!,  4  AS  N!(l),  4  AS  N!<2),  4  AS  N!<3),  4  AS  N!<4),  4  AS  N!<5),  4  AS  N!(6),  3  AS  N!(7),  4  AS  N$<8),  4  A 

S  N!(9) ,  4  AS  N!(10) ,  4  AS  N!<11),  8  AS  N!(12),  4  AS  N!(13),  4  AS  N!(14),  4  AS  N!(15) ,  4  AS  N!(16),  8  AS  N!(17) 

F1=PN!+ '.elements':  OPEN  F!  AS  112  LEN=94 

FIELDH2.6  AS  Flg2!,2  AS  Lt!,2  AS  Rt!,4  AS  E!<1),4  AS  E!(2),4  AS  E!<3),4  AS  E!(4),4  AS  E!(5),4  AS  E!(6),4  AS  E!(7) ,4  AS  E 

!(8) ,4  AS  E!(9) ,4  AS  E!<10),4  AS  E!(ll),4  AS  E!<12),8  AS  E!<13),4  AS  E!(14),4  AS  E!(15),4  AS  E!(16),4  AS  E!(17),4  AS  E!< 

18), 8  AS  E!( 19) 

F!=PN!+'. displ':  OPEN  F!  AS  113  LEN=24:  FIELDS, 8  AS  W,8  AS  V!,8  AS  Ac! 

Load. Global .Matrices:  BUTTON  1,2 

CALL  Retrieve .Matrix(n,n+1, KUO, PN!+' .K&F.c' ,n4) 

CALL  Retrieve.Matrix(n,n,MIK),PN!*,.M.c',n4) 

CALL  Retrieve .Matrix(n,n3,U0il<),PN!*’. ini tial' ,n4) 

FOR  i=l  TO  n  'start  'displ '  file  3  zero 
LSET  U!=MKD!(U0»(i,l)):LSET  V!=NKD!(U0ll(  i  ,2))  :LSET  Ac!=NKD!(U0l< i ,3)) 

PUT83, i :NEXT  i 
BUTTON  1,0 

'*«««mm**m»**«*s«m«**«**««*«««m  debug  only 
CALL  Displ ayMatrix(n,n*l,KI(), ’Stiffness*) 

CALL  Displ ayMatrix<n,n, MHO, 'Mass’) 

CALL  DisplayMatrix(n,n3,U0l(), 'Initial  Conditions') 

'Get  or  calculate  constants 
delta=l/2:alpha=l/4:  GOSUB  Get. deltaT. and. Time. Steps 

A0=l/(alpha»DeltaT*2):A2=l/(alpha*DeltaT):A3=l/(2»alpha)-l  '  calculate  constants 
A4=Del taT*(l-del ta) :A7=del ta*Del taT 

FOR  Counter=l  TO  NumSteps  '  begin  solution  loop 

Find.Dyn.Force.Nat:  BUTTON  2,2  '  status  report 
GOSUB  AssembleForceMat:  GOSUB  El enentMatrixAssembler 
BUTTON  2,1  '  status  report 

Find.Effective.Mat:  BUTTON  3,2  '  also  apply  BC 

FOR  i=l  TO  N:KII(i  ,n*l)=KI(i  ,n*  1  )*Fdl<  i ) :  Fdl<  i  )=0  'add  in  dyn  forces  and  in  it  Fdi  for  next  time  step 
IF  U1 1(1,1)01  THEN  FOR  j=l  TO  N:Kt( i ,N*1)=KI( i ,n*l)*MI( i ,j)e(AO*UOI( j ,l)*A2»U0t( j ,2)*A3»U0l(j ,3)) :KI< i ,  j)=KI< i , j)+AO* 
MKi  ,j):NEXT  j 
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IF  Ult(i ,1)31  THEN  GQSUB  Essential  .B.C:  FOR  j=l  TO  n :K#< i  f j)=-< i=j ) :NEXT  j :K»<  i  ,n+l)=Displ  'set  spec'd  displacement 
NEXT  i 

BUTTON  3,1  '  status  report 

Solve:  BUTTON  4,2  '  status  report 
60SUB  Guass:  BUTTON  4,1:  BUTTON  5,2 

FOR  i=l  TO  n  '  tind  V  and  A  vectors  and  store  displacements  in  U 
UW(i,l)=KI(i,nM) 

Ul«(  i  ,3)=A0*(Uli<  i ,  1  )-U0tt<  i  ,1))-A2*U0»<  i  ,2)-A3*U0<)<  i  ,3) 

Ul«(  i  ,2)=U0#<  i  ,2)+A4*U0i<  i  ,3)+A7*Ullt(  i  ,3) 

LSET  W=NKM(Ul#<i,l)):LSET  V4=NKI»<Ul«(i ,2)):LSET  Ac«=MKD$(Ul»< i , 3) > : j=  i ♦Counter*n :PUT«3, j  '  save  to  disk 
NEXT  i:  BUTTON  5,1  '  status  report 
' HIHHiHHHiHHUHiHiHHHHHHi  debug  Only 

CALL  Di spT aytiatr ix<n,n3,Ultt<), 'Displacement,  Velocity,  and  Acceleration1) 

UINDOU  2:  PRINT  USIN6  'Time  step  «»l  ot  m  .'jcounter.NuraSteps 
PRINT  USING  'T  3  »».»**"*  Time  step  =  «».«« . ;T,deltaT:  UINDOU  1 

'HHHHlHHHHHHIHHIIIHiHiHiii 

NextTimeStep:  BUTTON  4,2:  UINDOU  OUTPUT  2  '  status  report 
T=T+deltaT:  UINDOU  i 

FOR  i3l  TO  n:  FOR  j=l  TO  3:  U0*<i,j)3Ul»<i,j):  NEXT  j,i  '  intialize 
CALL  Retrieve .Natrix(n,n+l,KI<),PN$+*.K4F.c,ln4) 

BUTTON  4,1  'status  report 
NEXT  Counter 

CLOSE:  UINDOU  CLOSE  1:  UINDOU  CLOSE  2:  CHAIN  'DynFEP.menu'i  END 
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A  subset  ot  the  toll  owing  SUB-Programs  are  used  in  most  o4  the  DynFEP  programs: 


Sub-Programs  Below 


SUB  Retrieve.Matrix(r,c,AIO,FI,k)  STATIC 
IF  UBOllND(AI,lXr  OR  U80tND<AI,2Xc  THEN  PRINT  CHR*<?)*Fatal  error!1:  STOP 
RL=C*8:  OPBi  Fi  AS  Ik  LEN=RL:  FIELDHk ,RL  AS  AA$ 

FOR  i=l  TO  r:  GETIk , i :  FOR  j=l  TO  c 
BI=NIDI<#»,8*<  j-lXl  ,8) ;  Al<  i ,  j)=CW<BI) 

NEXT  j , i :  CLOSEIk 
BID  SUB 

SUB  Store  .Natrix(r,c,AIO, FI, k)  STATIC 
IF  UBOWD<AI,lXr  OR  UB0UND<AI,2Xc  THEN  PRINT  CHRI<7)*Fatal  error!1:  STOP 
RL=C*8:  OPEN  FI  AS  Ik  LEN=RL:  FIELOtk ,RL  AS  AAI 
FOR  i=l  TO  r:  81=' *:  FOR  j=l  TO  c 
BI=Bl4NKOKAX(i,j)) 

NEXT  j:  LSET  AAI=BI:  PUTXk.i:  NEXT  i:  CLOSE  Ik 
END  SUB 

SUB  D i sp 1  ay .Natr i x < Row , Col ,At ( 2) ,TI)  STATIC 
CALL  TEXTFONTU):  CALL  TEXTSI2E(9):  PRINT  Tl 
FOR  isl  TO  Row:  FOR  j=l  TO  Col 

PRINT  USING  M.H . ;AX<i,j)i 

NEXT  j:  PRINT:  NEXT  i:  PRINT 
INPUT  'Press  'RETURN'  to  continue1 ;al 
END  SUB 

SUB  Nat. times. Mat(rAlcB1cArB,AI(2)pBI<2),RI<2))  STATIC 
'IAI  »  183  =  [RI 

'rA  =  irows  in  IAI  cArB  =  Icols  in  CA3  and  Irows  in  EB3 

'cB  =  Icols  in  CB1  CR1  is  dimensioned  rA  X  cB 

FOR  i=l  TO  rA:  FOR  j=l  TO  cB:  FOR  k=l  TO  cArB 
Rl<  i ,  j)=RI(  i ,  jXAK  i  ,k)*BI(k ,  j) 

NEXT  k,j,i 
BID  SUB 

SUB  NatTrans.times.Hat<cAlcB,rArB,AI<2),BK2)lRI<2))  STATIC 
'IA  transpose!  •  CB3  =  (R3 

'cA  =  Icols  in  (A1  rArB  =  Irows  in  [A]  and  Irows  in  IBI 

'cB  =  Icols  in  (B]  CR3  is  dimensioned  cA  X  cB 

FOR  i=J  TO  cA:  FOR  j=l  TO  cB:  FOR  k=l  TO  rArB 
RKi  ,j)=RI(i  ,jXAI(k,i)*BI(k,j) 

NEXT  k,j,i 
BID  SUB 

SUB  Nat.p1us.Nat(r,c,Cll,AI(2),C2l,BI<2))  STATIC 
'Cl*tAl  ♦  C2«IB1  =  result  stored  in  tAl 

FOR  i=l  TO  r:  FOR  j*l  TO  c:  AI<i,j)=Cll*AI<i,jXC2MBKi,j>:  NEXT  j,i 
BID  SUB 

SUB  Invert.Hatr ix<n,AI<2))  STATIC 
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'Takes  IAJ  *  CAl'-l  =  [11  AND  changes  TO  Ill  *  £A3 *-l  =  IA1*-1  (based  on  Guass  elimination) 

'  IAJ‘-1  replaces  [A1 

DIN  ll<n,n):  FOR  i=l  TO  n:  Il(i,i)=l:  NEXT  i  'identity  matrix 

FOR  i=l  TO  n:  m#=AKi,i):  FOR  j=l  TO  n:  A»(i ,j)=AI(i ,j)/m»:  I«(i,j)=I«(i,j)/m»:  NEXT  j 

FOR  k=l  TO  n:  IF  KOi  THBt  n«=A«<k f i ) :  FOR  j=l  TO  n:  AI(k,j)=A»<k,j)-A*(i,j)*m»:  l*(k,j)=l»(k,j)-l*(i,j)*m»:  NEXT  j 
NEXT  k,i 

FOR  i=l  TO  n;  FOR  j=l  TO  n:  A*(i ,j)=I»(i ,j):  NEXT  j, i :  ERASE  Ik  '  store  inverse  in  A« 

END  SUB 

SUB  Determinants, AK2),Det)  STATIC 
'Uses  pivital  condensation  to  find  the  determinant  of  [All 
Nult=l:  Signal 
WHILE  n=>2 

i=l:  WHILE  A8(i,l)=0  AND  i<=n:  i= i + 1 :  WEND  'check  for  zero  in  first  column,  then  correct 
IF  i)n  THEN  Det=0 :  GOTO  Finished  '1st  col  has  all  zeros 

IF  i)l  THEN  SigiF-Sign:  FOR  j=l  TO  n:  SWAP  AI(i,j),Ak(i,j):  NEXT  j  'swap  rows  and  change  sign 
Mu l t=Nul t/AI( 1 , 1 ) ‘ ( n-2) 

FOR  i=2  TO  n:  FOR  j=2  TO  n:  Akti ,j)=A#(l,l)«A«(i ,j)-AI(l,j)*A«(i ,1):  NEXT  j,i 
FOR  i=l  TO  n-ls  FOR  j=l  TO  n-ls  A»<i,j)=AI<iel,j+l):  NEXT  j,i 
n=n-l 
WEND 

Det=Sign*Mu1tsAI(l,l) 

Finished: 
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